x = 3
y = 1.2
z = "Hi!"
"Hi!"
typeof(x)
Int64
typeof(y)
Float64
typeof(z)
String
x + 1
4
# \mu + tab
μ = 0.0
σ = 2.0
xᵢ = 0
σ² = σ^2
4.0
x = 1.0
y = 2.0
if x < y
println("x is less than y.")
elseif x == y
println("x is equal to y.")
else
println("x is greter than y.")
end
x is less than y.
x = -3
# x < 0が成り立てば:の手前をyに代入する。成り立たねければ:の後をyに代入する
y = x < 0 ? "T" : "F"
println(y)
T
a = 0
b = 0
true && (a = 1) # (A && B)AがtrueならばBを実行
true || (b = 1) # (A || B)AがfalseならばBを実行
println("a = $(a), b = $(b)")
a = 1, b = 0
for i in 1:10
println(i)
end
1 2 3 4 5 6 7 8 9 10
myinv(x) = 1 / x
myinv (generic function with 1 method)
myinv(3)
0.3333333333333333
function mymean(x, y)
return (x + y) / 2
end
mymean (generic function with 1 method)
mymean(1.0, 2.0)
1.5
# 縦ベクトル
a = [1, 2, 3]
3-element Vector{Int64}:
1
2
3
# 横ベクトル
b = [1 2 3]
1×3 Matrix{Int64}:
1 2 3
c = Array{Float64}(undef, 3)
3-element Vector{Float64}:
2.2439480934e-314
2.263943333e-314
2.2426095865e-314
d = zeros(3)
3-element Vector{Float64}:
0.0
0.0
0.0
e = ones(3)
3-element Vector{Float64}:
1.0
1.0
1.0
# 0~1一様
f = rand(3)
3-element Vector{Float64}:
0.31385851194472925
0.9294287946579434
0.5799688477670801
# 標準正規分布
g = randn(3)
3-element Vector{Float64}:
-2.1397517175103573
-0.3955230987292779
0.31163264486238834
A = [1 2 3 4;
5 6 7 8]
2×4 Matrix{Int64}:
1 2 3 4
5 6 7 8
B = ones(2,4)
2×4 Matrix{Float64}:
1.0 1.0 1.0 1.0
1.0 1.0 1.0 1.0
size(A)
(2, 4)
length(A)
8
A[2,1]
5
A[:,1]
2-element Vector{Int64}:
1
5
A[:,1:3]
2×3 Matrix{Int64}:
1 2 3
5 6 7
[2*i for i in 1:5]
5-element Vector{Int64}:
2
4
6
8
10
[i*j for i in 1:3, j in 1:4]
3×4 Matrix{Int64}:
1 2 3 4
2 4 6 8
3 6 9 12
params1 = (1, 2, 3)
(1, 2, 3)
f1(a, b, c) = a + b + c
f1(params1...)
6
params2 = (1, 2, 3, 4)
(1, 2, 3, 4)
f2(a, b, c, d) = a + b + c + d
f2(params2...)
10
a = [1, 2, 3]
a + 1
MethodError: no method matching +(::Vector{Int64}, ::Int64)
For element-wise addition, use broadcasting with dot syntax: array .+ scalar
Closest candidates are:
+(::Any, ::Any, ::Any, ::Any...) at operators.jl:560
+(::T, ::T) where T<:Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8} at int.jl:87
+(::Base.TwicePrecision, ::Number) at twiceprecision.jl:267
...
Stacktrace:
[1] top-level scope
@ In[35]:2
[2] eval
@ ./boot.jl:360 [inlined]
[3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1116
a .+ 1
3-element Vector{Int64}:
2
3
4
# ブロードキャスト
function add_one(x)
x + 1
end
add_one(1.0)
2.0
add_one.([1,2,3])
3-element Vector{Int64}:
2
3
4
x -> x + 1
#5 (generic function with 1 method)
map(x -> x+1, [1,2,3])
3-element Vector{Int64}:
2
3
4
function test(maxiter)
a = []
for i in 1:maxiter
push!(a, randn())
end
sum(a)
end
@time test(100000)
0.008677 seconds (201.50 k allocations: 5.155 MiB, 36.59% compilation time)
165.426717533415
mean([1,2,3,4,5])
UndefVarError: mean not defined Stacktrace: [1] top-level scope @ In[42]:1 [2] eval @ ./boot.jl:360 [inlined] [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String) @ Base ./loading.jl:1116
"""
よく使うパッケージ
Distirbutions.jl
ForwordDiff.jl
IJulia.jl
PyPlot.jl
"""
using Pkg
Pkg.add("Statistics")
cannot document the following expression: using Pkg Stacktrace: [1] error(::String, ::String) @ Base ./error.jl:42 [2] top-level scope @ In[43]:1 [3] eval @ ./boot.jl:360 [inlined] [4] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String) @ Base ./loading.jl:1116
using Statistics
mean([1,2,3,4,5])
3.0
Pkg.add("PyPlot")
UndefVarError: Pkg not defined Stacktrace: [1] top-level scope @ In[45]:1 [2] eval @ ./boot.jl:360 [inlined] [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String) @ Base ./loading.jl:1116
using PyPlot
# 1つ目のデータ
Y1 = [1, 7, 11, 13, 15, 16]
# 2つ目のデータ
Y2 = [15, 3, 13, 2, 7, 1]
# グラフを作成
fig, ax = subplots()
# 折れ線グラフを描画
ax.plot(Y1)
ax.plot(Y2)
# 軸の名称を設定
ax.set_xlabel("index")
ax.set_ylabel("Y")
# グラフのタイトルを設定
ax.set_title("my graph")
PyObject Text(0.5, 1.0, 'my graph')
fig, axes = subplots(1, 2, figsize=(10,3))
axes[1].plot(Y1)
axes[2].plot(Y2)
1-element Vector{PyCall.PyObject}:
PyObject <matplotlib.lines.Line2D object at 0x7fe074f915e0>
# ２次関数の定義
quadraticf(x) = x^2
# -3から3まで、10個の等間隔の点列を生成
xs = range(-3, 3, length=10)
fig, ax = subplots()
# 関数の描画
ax.plot(xs, quadraticf.(xs), "-")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title("y=x^2")
PyObject Text(0.5, 1.0, 'y=x^2')
# -3から3まで、10個の等間隔の点列を生成
xs = range(-3, 3, length=100)
fig, ax = subplots()
# 関数の描画
ax.plot(xs, quadraticf.(xs), "-")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title("y=x^2")
PyObject Text(0.5, 1.0, 'y=x^2')
r = 1.0
fx(θ) = r*(θ-sin(θ))
fy(θ) = r*(1-cos(θ))
θs = range(-3pi, 3pi, length=100)
fig, ax = subplots()
ax.plot(fx.(θs), fy.(θs))
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title("cycloid")
PyObject Text(0.5, 1.0, 'cycloid')
# シグモイド関数を定義
sig(x) = 1 / (1 + exp(-x))
xs = range(-5, 5, length=100)
fig, ax = subplots()
ax.plot(xs, sig.(xs))
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title("sigmoid")
PyObject Text(0.5, 1.0, 'sigmoid')
# 可視化に使う乱数を生成
X = randn(1000)
Y = randn(1000)
fig, ax =subplots()
# ax.plot(X, Y, "o")などでも可
ax.scatter(X, Y)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title("scatter")
PyObject Text(0.5, 1.0, 'scatter')
X = randn(1000)
fig, ax = subplots()
ax.hist(X)
([2.0, 22.0, 47.0, 145.0, 212.0, 246.0, 186.0, 86.0, 39.0, 15.0], [-3.3668178997339817, -2.7338823046905456, -2.1009467096471095, -1.468011114603673, -0.8350755195602368, -0.20213992451680074, 0.4307956705266358, 1.0637312655700715, 1.696666860613508, 2.3296024556569446, 2.96253805070038], (PyObject <matplotlib.patches.Rectangle object at 0x7fe05f453af0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f444b80>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f453fa0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f4641f0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f464400>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f464610>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f464820>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f464a30>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f464c40>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f464e50>))
X = randn(1000)
fig, ax = subplots()
ax.hist(X, bins=50)
([3.0, 0.0, 2.0, 1.0, 0.0, 3.0, 4.0, 4.0, 5.0, 7.0 … 3.0, 5.0, 2.0, 1.0, 1.0, 3.0, 0.0, 0.0, 0.0, 1.0], [-3.3940216424986653, -3.2567722674011526, -3.1195228923036398, -2.9822735172061265, -2.8450241421086138, -2.707774767011101, -2.5705253919135878, -2.433276016816075, -2.296026641718562, -2.1587772666210494 … 2.233202736499363, 2.370452111596876, 2.5077014866943887, 2.6449508617919024, 2.782200236889415, 2.919449611986928, 3.0566989870844408, 3.1939483621819535, 3.3311977372794663, 3.46844711237698], (PyObject <matplotlib.patches.Rectangle object at 0x7fe05f60deb0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f4dff40>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f61c3a0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f61c5b0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f61c7c0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f61c9d0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f61cbe0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f61cdf0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f61cfa0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f62d250>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f62d460>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f62d670>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f62d880>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f62da90>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f62dca0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f62deb0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f63b100>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f63b310>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f63b520>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f63b730>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f63b940>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f63bb50>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f63bd60>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f63bf70>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f64b1c0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f64b3d0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f64b5e0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f64b7f0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f64ba00>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f64bc10>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f64be20>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f64bfd0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f65b280>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f65b490>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f65b6a0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f65b8b0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f65bac0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f65bcd0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f65bee0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f668130>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f668340>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f668550>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f668760>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f668970>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f668b80>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f668d90>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f668fa0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f6791f0>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f679400>, PyObject <matplotlib.patches.Rectangle object at 0x7fe05f679610>))
# 2変数関数を定義
fz(x, y) = exp(-(2x^2+y^2+x*y))
# x軸とy軸の可視化範囲を定義
xs = range(-1,1, length=100)
ys = range(-2,2, length=100)
fig, axes = subplots(1,2, figsize=(8,4))
# 1つ目のsubplotには数値をつける
cs = axes[1].contour(xs, ys, fz.(xs', ys))
axes[1].clabel(cs, inline=true)
# 2つ目のsubplotにはカラーバーをつける
cs = axes[2].contourf(xs, ys, fz.(xs', ys))
fig.colorbar(cs)
# subplots間の余白の大きさを自動調整
tight_layout()
a = [1, 2, 3]
2*a
3-element Vector{Int64}:
2
4
6
b = [4, 5, 6]
a + b
3-element Vector{Int64}:
5
7
9
a .* b
3-element Vector{Int64}:
4
10
18
#内積
sum(a .* b)
32
a' * b
32
A = [1 2 3;
4 5 6]
B = [10 20 30;
40 50 60]
A + B
2×3 Matrix{Int64}:
11 22 33
44 55 66
A = [1 2;
3 4;
5 6]
B = [10 20 30 40;
50 60 70 80]
C = A * B
3×4 Matrix{Int64}:
110 140 170 200
230 300 370 440
350 460 570 680
M = size(A,1) #Aの行数
N = size(B, 2) #Bの列数
# M * Nの行列を作成
C = [sum(A[i,:] .* B[:,j]) for i in 1:M, j in 1:N]
3×4 Matrix{Int64}:
110 140 170 200
230 300 370 440
350 460 570 680
B * A
DimensionMismatch("matrix A has dimensions (2,4), matrix B has dimensions (3,2)")
Stacktrace:
[1] _generic_matmatmul!(C::Matrix{Int64}, tA::Char, tB::Char, A::Matrix{Int64}, B::Matrix{Int64}, _add::LinearAlgebra.MulAddMul{true, true, Bool, Bool})
@ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/matmul.jl:814
[2] generic_matmatmul!(C::Matrix{Int64}, tA::Char, tB::Char, A::Matrix{Int64}, B::Matrix{Int64}, _add::LinearAlgebra.MulAddMul{true, true, Bool, Bool})
@ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/matmul.jl:802
[3] mul!
@ /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/matmul.jl:302 [inlined]
[4] mul!
@ /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/matmul.jl:275 [inlined]
[5] *(A::Matrix{Int64}, B::Matrix{Int64})
@ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/matmul.jl:153
[6] top-level scope
@ In[65]:1
[7] eval
@ ./boot.jl:360 [inlined]
[8] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1116
A = [1 2;
3 4;
5 6]
I = [1 0 0;
0 1 0;
0 0 1]
I * A
3×2 Matrix{Int64}:
1 2
3 4
5 6
A = [1 2 3;
4 5 6]
A'
3×2 adjoint(::Matrix{Int64}) with eltype Int64:
1 4
2 5
3 6
a = [1, 2, 3]
b = [5, 7]
a * b'
3×2 Matrix{Int64}:
5 7
10 14
15 21
f2(x, y) = 2*x + y
f2.(a, b')
3×2 Matrix{Int64}:
7 9
9 11
11 13
A = [1 2;
3 4]
B = inv(A)
2×2 Matrix{Float64}:
-2.0 1.0
1.5 -0.5
A * B
2×2 Matrix{Float64}:
1.0 0.0
8.88178e-16 1.0
B * A
2×2 Matrix{Float64}:
1.0 0.0
2.22045e-16 1.0
println(B)
[-1.9999999999999996 0.9999999999999998; 1.4999999999999998 -0.4999999999999999]
A = Rational{Int}[1 2;
3 4]
B = inv(A)
println(A * B)
println(B * A)
Rational{Int64}[1//1 0//1; 0//1 1//1]
Rational{Int64}[1//1 0//1; 0//1 1//1]
A = Rational{Int}[1 2;
3 4]
sol = inv(A) * [-1, 1]
2-element Vector{Rational{Int64}}:
3//1
-2//1
using Statistics
X = rand(5)
5-element Vector{Float64}:
0.19604681062925633
0.15039929085606318
0.6494316784702181
0.5041200110315147
0.2762255348491489
Y = rand(2, 5)
2×5 Matrix{Float64}:
0.880842 0.610478 0.409949 0.1643 0.516474
0.488594 0.386113 0.834355 0.291409 0.122172
println(sum(X))
println(mean(X))
1.7762233258362012 0.35524466516724024
println(sum(Y))
println(sum(Y, dims=1))
println(sum(Y, dims=2))
4.704685348250634 [1.3694354496341408 0.996590834143914 1.2443035609926136 0.4557090946317597 0.6386464088482056] [2.582042658726274; 2.1226426895243593]
println(std(X))
println(std(X)^2)
println(var(X))
0.21345929039558376 0.045564868656186155 0.045564868656186155
cov(Y,dims=1)
5×5 Matrix{Float64}:
0.0769292 0.0440033 -0.0832362 -0.024929 0.0773321
0.0440033 0.0251697 -0.0476108 -0.0142593 0.0442337
-0.0832362 -0.0476108 0.0900602 0.0269728 -0.0836721
-0.024929 -0.0142593 0.0269728 0.00807828 -0.0250596
0.0773321 0.0442337 -0.0836721 -0.0250596 0.077737
cov(Y,dims=2)
2×2 Matrix{Float64}:
0.0692436 0.00573912
0.00573912 0.0706695
Pkg.add("Distributions")
UndefVarError: Pkg not defined Stacktrace: [1] top-level scope @ In[84]:1 [2] eval @ ./boot.jl:360 [inlined] [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String) @ Base ./loading.jl:1116
using Distributions
μ = 1.5
σ = 2.0
Z = rand(Normal(μ, σ), 10000)
Z
10000-element Vector{Float64}:
-2.878438829554005
1.9253341325276128
0.7912537284479731
-0.28884382408978526
7.4798218320469
1.4923605551865857
3.0248292221970154
1.039516574189463
5.394882281609656
0.14356031974931138
0.4469234226249874
1.8896840550067933
3.5436087455003693
⋮
-0.14534552137994927
-1.7913541860662607
2.4454416739331584
0.007851113359583328
-0.6249474431000608
2.12050537813308
-1.0890719963791056
-1.279649551121989
0.9741969853930287
2.583001380041292
1.5100434409297825
0.34536356361712217
println(mean(Z))
println(std(Z))
1.524948267501824 1.9954566091407018
Pkg.add("PyPlot")
UndefVarError: Pkg not defined Stacktrace: [1] top-level scope @ In[87]:1 [2] eval @ ./boot.jl:360 [inlined] [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String) @ Base ./loading.jl:1116
using PyPlot
# quadraticf(x)を２次関数として定義
quadraticf(x) = -(x + 1)*(x-1)
# hを微小な量として定義（10のマイナス10乗）
h = 1.0e-10
# 導関数quadraticf'の近似式
quadraticff(a) = (quadraticf(a+h) - quadraticf(a)) / h
# 関数の可視化範囲
xs = range(-1, 1, length=100)
fig, axes = subplots(2,1, figsize=(4,6))
# 関数のプロット
axes[1].plot(xs, quadraticf.(xs), "b")
axes[1].grid()
axes[1].set_xlabel("x"), axes[1].set_ylabel("y")
axes[1].set_title("function f(x)")
# 導関数のプロット
axes[2].plot(xs, quadraticff.(xs), "r")
axes[2].grid()
axes[2].set_xlabel("x"), axes[2].set_ylabel("y")
axes[2].set_title("derivative f'(x)")
tight_layout()
# グラフを可視化する解像度
L = 10
# f2(x)を可視化する範囲
xs1 = range(-1, 1, length=L)
xs2 = range(-1, 1, length=L)
# 2変数関数の定義
f2(x) = -(x .+ 1)' * (x .- 1)
# 2変数関数の勾配
∇f2(x) = -2x
fig, axes = subplots(1,2, figsize=(8,4))
# 関数の等高線の可視化
cs = axes[1].contour(xs1, xs2, [f2([x1,x2]) for x1 in xs1, x2 in xs2]')
axes[1].clabel(cs, inline=true)
axes[1].set_xlabel("x1"), axes[1].set_ylabel("x2")
axes[1].set_title("f2(x)")
# 勾配ベクトルの計算と可視化
vec1 = [∇f2([x1,x2])[1] for x1 in xs1, x2 in xs2]
vec2 = [∇f2([x1,x2])[2] for x1 in xs1, x2 in xs2]
axes[2].quiver(repeat(xs1, 1, L), repeat(xs2, L, 1), vec1, vec2)
axes[2].set_xlabel("x1"), axes[2].set_ylabel("x2")
axes[2].set_title("∇f2(x)")
tight_layout()
Pkg.add("ForwardDiff")
UndefVarError: Pkg not defined Stacktrace: [1] top-level scope @ In[90]:1 [2] eval @ ./boot.jl:360 [inlined] [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String) @ Base ./loading.jl:1116
using ForwardDiff
# 2変数関数の定義
f3(x) = -(x + 1) * (x - 1)
# 自動微分によって導関数f'(x)を求める
f33(x) = ForwardDiff.derivative(f3, x)
# 関数の可視化範囲
xs = range(-1, 1, length=100)
fig, axes = subplots(2,1, figsize=(4,6))
# 関数のプロット
axes[1].plot(xs, f3.(xs), "b")
axes[1].grid()
axes[1].set_xlabel("x"), axes[1].set_ylabel("y")
axes[1].set_title("functin f(x)")
# 導関数のプロット
axes[2].plot(xs, f33.(xs), "r")
axes[2].grid()
axes[2].set_xlabel("x"), axes[1].set_ylabel("y")
axes[2].set_title("derivative f'(x)")
tight_layout()
fig, ax = subplots()
xs = range(0, 2pi*3, length=100)
# sin(x)をプロット
ax.plot(xs, sin.(xs), color="b", label="sin(x)")
# 導関数をプロット
ax.plot(xs, map(x -> ForwardDiff.derivative(sin, x), xs), color="r", label="sin'(x)")
ax.grid()
ax.set_xlabel("x"), ax.set_ylabel("y")
ax.legend()
PyObject <matplotlib.legend.Legend object at 0x7fe0608ff310>
# シグモイド関数を定義
sig(x) = 1/(1 + exp(-x))
xs = range(-5, 5, length=100)
fig, axes = subplots(2, 1, figsize=(4,6))
# シグモイド関数をプロット
axes[1].plot(xs, sig.(xs), color="b")
axes[1].set_xlabel("x"), axes[1].set_ylabel("y")
axes[1].set_title("sig(x)")
axes[1].grid()
# 導関数をプロット
axes[2].plot(xs, map(x -> ForwardDiff.derivative(sig, x), xs), color="r")
axes[2].set_xlabel("x"), axes[2].set_ylabel("y")
axes[2].set_title("sig'(x)")
axes[2].grid()
# 最大値を探したい関数
x_opt = 0.50
quadraticf(x) = -2(x - x_opt)^2
fig, ax = subplots()
xs = range(-3, 3, length=100)
ax.plot(xs, quadraticf.(xs), label="function")
ax.plot(x_opt, quadraticf(x_opt), "rx", label="optimal")
ax.set_xlabel("x"), ax.set_ylabel("y")
ax.grid()
ax.legend()
PyObject <matplotlib.legend.Legend object at 0x7fe060e2c3a0>
# 1変数関数の最適化
function gradient_method_1dim(f, x_init, η, maxiter)
# 最適化過程のパラメータを格納する配列
x_seq = Array{typeof(x_init), 1}(undef, maxiter)
# 勾配
ff(x) = ForwardDiff.derivative(f, x)
# 初期値
x_seq[1] = x_init
# メインの最適解ループ
for i in 2:maxiter
x_seq[i] = x_seq[i-1] + η*ff(x_seq[i-1])
end
x_seq
end
gradient_method_1dim (generic function with 1 method)
# 探索範囲の初期値
x_init = -2.5
# 探索の繰り返し数
maxiter = 20
# ステップサイズ
η = 0.1
# 最適化計算を実行
x_seq = gradient_method_1dim(quadraticf, x_init, η, maxiter)
quadraticf_seq = quadraticf.(x_seq)
20-element Vector{Float64}:
-18.0
-6.479999999999999
-2.3327999999999993
-0.8398079999999998
-0.30233087999999986
-0.10883911679999994
-0.03918208204799999
-0.014105549537279988
-0.0050779978334207915
-0.0018280792200314835
-0.0006581085192113332
-0.0002369190669160795
-8.529086408978862e-5
-3.07047110723239e-5
-1.1053695986036396e-5
-3.979330554973165e-6
-1.4325589997903394e-6
-5.157212399245448e-7
-1.8565964637284966e-7
-6.683747269421776e-8
# 目的関数の値をステップごとにプロット
fig, ax = subplots(figsize=(8,3))
ax.plot(quadraticf_seq)
ax.set_xlabel("iteration"), ax.set_ylabel("f")
ax.grid()
fig,axes = subplots(1,2,figsize=(8,3))
# 関数のプロット
axes[1].plot(xs, quadraticf.(xs), label="function")
# 探索の過程
axes[1].plot(x_seq, quadraticf.(x_seq), "bx", label="sequence")
# 最適値
axes[1].plot(x_opt, quadraticf(x_opt), "rx", label="optimal")
axes[1].set_xlabel("x"), axes[1].set_ylabel("y")
axes[1].grid()
axes[1].legend()
# 探索の過程をステップごとにプロット
axes[2].plot(1:maxiter, x_seq, "bx-", label="sequence")
axes[2].hlines(x_opt, 0, maxiter, ls="--", label="x_opt")
axes[2].set_xlim([0, maxiter])
axes[2].set_xlabel("iteration"), axes[2].set_ylabel("x")
axes[2].grid()
axes[2].legend()
tight_layout()
# ２次関数を定義
x_opt = [0.50, 0.25]
f2(x) = -sqrt(0.05 + (x[1] - x_opt[1])^2) - (x[2] - x_opt[2])^2
# 関数を等高線として可視化
L = 100
xs1 = range(-1, 1, length=L)
xs2 = range(-1, 1, length=L)
fig, ax = subplots()
ax.contour(xs1, xs2, [f2(x1, x2) for x1 in xs1, x2 in xs2]')
ax.plot(x_opt[1], x_opt[2], color="r", marker="x", label="optimal")
ax.set_xlabel("x1"), ax.set_ylabel("x2")
ax.grid()
ax.legend()
PyObject <matplotlib.legend.Legend object at 0x7fe060f2dd00>
# 多変数関数のための勾配法
function gradient_method(f, x_init, η, maxiter)
# 最適化過程のパラメータを格納する配列
x_seq = Array{typeof(x_init[1]), 2}(undef, length(x_init), maxiter)
# 勾配
∇f(x) = ForwardDiff.gradient(f, x)
# 初期値
x_seq[:, 1] = x_init
# メインの最適化ループ
for i in 2:maxiter
x_seq[:, i] = x_seq[:, i-1] + η*∇f(x_seq[:, i-1])
end
x_seq
end
# パラメータの設定
x_init = [-0.75, -0,75]
maxiter = 20
η=0.1
# 最適化の実行
x_seq = gradient_method(f2, x_init, η, maxiter)
f_seq = [f2(x_seq[:,i]) for i in 1:maxiter]
20-element Vector{Float64}:
-1.3323425099200294
-1.2130713587869222
-1.10246741292214
-0.9977733058718217
-0.8973440357566991
-0.8003014554196972
-0.7063454734180662
-0.6156878830524037
-0.5291020764458453
-0.448101220575972
-0.3752360053913261
-0.3143005037996588
-0.26957995148648584
-0.24287708480788592
-0.23050261870859573
-0.2258762089573429
-0.22433899880772853
-0.22384834870321288
-0.22369132580962012
-0.22363942129961714
# 目的関数の値をステップごとにプロット
fig, ax = subplots(figsize=(8,3))
ax.plot(f_seq)
ax.set_xlabel("iteration"), ax.set_ylabel("f")
ax.grid()
fig, axes = subplots(1, 2, figsize=(8, 3))
# 等高線図で関数を可視化
axes[1].contour(xs1, xs2, [f2(x1, x2) for x1 in xs1, x2 in xs2]')
# 最適化の過程
axes[1].plot(x_seq[1, :], x_seq[2,:], ls="--", marker="x", label="sequences")
axes[1].plot(x_opt[1], x_opt[2], color="r", marker="x", label="optimal")
axes[1].set_xlabel("x1"), axes[1].set_ylabel("x2")
axes[1].grid()
axes[1].legend()
# ステップごとの最適化の過程
axes[2].plot(1:maxiter, x_seq[1,:], color="b", marker="o", label="x[1]")
axes[2].plot(1:maxiter, x_seq[2,:], color="r", marker="^", label="x[2]")
axes[2].hlines(x_opt[1], 0, maxiter, color="b", ls="--", label="x_opt[1]")
axes[2].hlines(x_opt[2], 0, maxiter, color="r", ls="-", label="x_opt[2]")
axes[2].set_xlabel("iteration")
axes[2].set_ylabel("x1, x2")
axes[2].grid()
axes[2].legend()
tight_layout()
# 目的関数の定義
f_complex(x) = 0.3*cos(3pi*x) - x^2
# 最適化のパラメータ設定
maxiter = 20
η = 0.01
# 初期値a
x_init_a = -0.8
x_seq_a = gradient_method_1dim(f_complex, x_init_a, η, maxiter)
f_seq_a = f_complex.(x_seq_a)
# 初期値b
x_init_b = 0.25
x_seq_b = gradient_method_1dim(f_complex, x_init_b, η, maxiter)
f_seq_b = f_complex.(x_seq_b)
# 目的関数の値をステップごとにプロット
fig, axes = subplots(2, 1, figsize=(8,4))
axes[1].plot(f_seq_a)
axes[1].set_xlabel("iteration"), axes[1].set_ylabel("f")
axes[1].set_title("f(x) (a)")
axes[1].grid()
axes[2].plot(f_seq_b)
axes[2].set_xlabel("iteration"), axes[2].set_ylabel("f")
axes[2].set_title("f(x) (b)")
axes[2].grid()
tight_layout()
# 関数を可視化する範囲
xs = range(-1, 1, length=100)
# 最適化の過程
fig, axes = subplots(1, 2, figsize=(12,4))
axes[1].plot(xs, f_complex.(xs), label="f_complex")
axes[1].plot(x_seq_a, f_complex.(x_seq_a), color="b", marker="x", label="x sequences")
axes[1].set_xlabel("x"), axes[1].set_ylabel("y")
axes[1].grid()
axes[1].set_title("initial value a")
axes[1].legend()
axes[2].plot(xs, f_complex.(xs), label="f_complex")
axes[2].plot(x_seq_b, f_complex.(x_seq_b), color="b", marker="x", label="x sequences")
axes[2].set_xlabel("x"), axes[2].set_ylabel("y")
axes[2].grid()
axes[2].set_title("initial value b")
axes[2].legend()
tight_layout()
# 学習用の入力値集合
X_obs = [1, 2, 4]
# 学習用の出力値集合
Y_obs = [10, 35, 72]
# データの可視化
fig, ax = subplots()
ax.scatter(X_obs, Y_obs)
ax.set_xlabel("x"), ax.set_ylabel("y")
ax.set_title("data")
ax.grid()
# 適当な傾きパラメータおよび切片パラメータを設定
w = [12.0, 10.0]
# 予測に使う関数
ff(x) = w[1]*x + w[2]
# 関数を可視化する範囲
xs = range(0, 5, length=100)
# データと予測関数の可視化
fig, ax = subplots()
ax.plot(xs, ff.(xs), label="function")
ax.scatter(X_obs, Y_obs)
ax.set_xlabel("x"), ax.set_ylabel("y")
ax.grid()
ax.legend()
PyObject <matplotlib.legend.Legend object at 0x7fe062cb7a00>
# 誤差関数の定義
E(w) = sum([(Y_obs[n] - (w[1]*X_obs[n] + w[2]))^2 for n in 1:length(X_obs)])
E (generic function with 1 method)
# 最適化するパラメータの初期値
w_init = [0.0, 0.0]
# 最適化計算の回数
maxiter = 500
#学習率
η = 0.01
# 最適化の実行、最大化アルゴリズムなので、-Eを目的関数にする
F(w) = -E(w)
w_seq = gradient_method(F, w_init, η, maxiter)
f_seq = [F(w_seq[:,i]) for i in 1:maxiter]
500-element Vector{Float64}:
-6509.0
-1939.5819999999997
-644.3813987199998
-275.7910161687038
-169.46811483893487
-137.41138630873562
-126.41782512660872
-121.4305759666068
-118.18766202638302
-115.48193910495918
-112.97076516030776
-110.556165388921
-108.20942661889302
⋮
-12.072442904302092
-12.072418776699926
-12.072395223013139
-12.072372229590176
-12.072349783104178
-12.072327870545308
-12.072306479213296
-12.072285596709893
-12.072265210931684
-12.072245310063245
-12.072225882570269
-12.072206917192597
w1, w2 = w_seq[:, end]
println("w_1 = $(w1), w_2 = $(w2)")
w_1 = 20.345436819385633, w_2 = -8.465882327774901
# 目的関数の値をステップごとにプロット
fig, ax = subplots(figsize=(8,4))
ax.plot(f_seq)
ax.set_xlabel("iteration"), ax.set_ylabel("f")
ax.grid()
# 予測に使う関数
ff(x) = w1*x + w2
# 関数を可視化する範囲
xs = range(0, 5, length=100)
# データと予測関数をプロット
fig, ax = subplots()
ax.plot(xs, ff.(xs), label="f(x)")
ax.scatter(X_obs, Y_obs, label="data")
ax.set_xlabel("x"), ax.set_ylabel("y")
ax.grid()
ax.legend()
PyObject <matplotlib.legend.Legend object at 0x7fe0615cef10>
function linear_fit(Y, X)
N = length(Y)
w1 = sum((Y .- mean(Y)) .* X) / sum((X .- mean(X)).*X)
w2 = mean(Y) - w1*mean(X)
w1, w2
end
w1, w2 = linear_fit(Y_obs, X_obs)
println("w1 = $(w1), w2 = $(w2)")
w1 = 20.35714285714286, w2 = -8.500000000000014
function approx_integration(x_range, f)
# 幅
Δ = x_range[2] - x_range[1]
# 近似された面積と幅を返す
sum([f(x) * Δ for x in x_range]), Δ
end
approx_integration (generic function with 1 method)
# 面積を近似したい関数
ff(x) = exp(-x^2)
# 等間隔の配列を用意
#x_range = range(-2, 2, length=10)
x_range = range(-5, 5, length=100)
# 積分近似の実行
approx, Δ = approx_integration(x_range, ff)
# 近似値（apporox）と厳密値（exact）の比較
println("approx = $(approx)")
println("exact = $(sqrt(pi))")
# 関数と近似結果の可視化
fig, ax = subplots(figsize=(8,4))
xs = range(-5, 5, length=1000)
ax.plot(xs, ff.(xs), "r--")
for x in x_range
ax.fill_between([x - Δ/2, x + Δ/2], [ff(x), ff(x)], [0, 0], facecolor="b", edgecolor="k", alpha=0.5)
end
# ax.set_xlabel("x"), ax.set_ylabel("y")
# ax.grid()
approx = 1.7724538509039651 exact = 1.7724538509055159
# 積分対象の２変数関数
D = 2
A = [0.5 0.3
0.3 1.0]
f2(x) = exp(-0.5*x'A*x)
# 20x20の区画に分割
L = 20
xs1 = range(-5, 5, length=L)
xs2 = range(-5, 5, length=L)
fig, axes = subplots(1, 2, figsize=(8,4))
# 等高線図で可視化
cs = axes[1].contour(repeat(xs1, 1, L), repeat(xs2', L, 1), [f2([x1, x2]) for x1 in xs1, x2 in xs2]')
axes[1].clabel(cs)
axes[1].grid()
axes[1].set_xlabel("x1"), axes[1].set_ylabel("x2")
# カラー表示
cs = axes[2].imshow([f2([x1, x2]) for x1 in xs1, x2 in xs2], origin="lower")
fig.colorbar(cs)
axes[2].set_xlabel("x1"), axes[2].set_ylabel("x2")
tight_layout()
function approx_integration_2dim(x_range, f)
Δ = x_range[2] - x_range[1]
sum([f([x1, x2]) * Δ^2 for x1 in x_range, x2 in x_range]), Δ
end
approx_integration_2dim (generic function with 1 method)
# det関数を使うためにLinearAlgebraパッケージを使用
Pkg.add("LinearAlgebra")
using LinearAlgebra
# 20x20の区間に分割
L = 20
x_range = range(-5, 5, length=L)
approx, Δ = approx_integration_2dim(x_range, f2)
println("approx = $(approx)")
println("exact = $(sqrt((2*pi)^D/det(A)))")
UndefVarError: Pkg not defined Stacktrace: [1] top-level scope @ In[114]:2 [2] eval @ ./boot.jl:360 [inlined] [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String) @ Base ./loading.jl:1116
# 1000x1000の区間に分割
L = 1000
x_range = range(-100, 100, length=L)
approx, Δ = approx_integration_2dim(x_range, f2)
println("approx = $(approx)")
println("exact = $(sqrt((2*pi)^D/det(A)))")
approx = 9.812686860654521
UndefVarError: det not defined Stacktrace: [1] top-level scope @ In[115]:7 [2] eval @ ./boot.jl:360 [inlined] [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String) @ Base ./loading.jl:1116
using Distributions
# パラメータが0.5のベルヌーイ分布を定義
bern = Bernoulli(0.5)
# 乱数を10個生成
X = rand(bern, 10)
10-element Vector{Bool}:
1
1
1
0
1
1
0
0
1
1
#パラメータを変更
bern = Bernoulli(0.9)
X = rand(bern, 10)
10-element Vector{Bool}:
1
1
1
1
1
1
0
0
0
1
bag(x::Bool) = x == 1 ? "A" : "B"
ball(y::Bool) = y == 1 ? "赤" : "白"
X = bag.(rand(bern, 10))
10-element Vector{String}:
"A"
"A"
"A"
"A"
"A"
"B"
"A"
"A"
"A"
"A"
function sample()
# 袋 の選択はそれぞれ1/2の確率
x = bag.(rand(Bernoulli(1//2)))
# 袋がAであれば、赤玉が出る確率は1/5
# 袋がBであれば、赤玉が出る確率は3/5
μ = x=="A" ? 1//5 : 3//5
# 玉の抽出
y = ball.(rand(Bernoulli(μ)))
x, μ, y
end
sample (generic function with 1 method)
for _ in 1:10
x, μ, y = sample()
println("袋: $(x), 玉: $(y)")
end
袋: A, 玉: 白 袋: B, 玉: 赤 袋: A, 玉: 赤 袋: A, 玉: 赤 袋: B, 玉: 赤 袋: A, 玉: 赤 袋: B, 玉: 白 袋: B, 玉: 赤 袋: A, 玉: 白 袋: B, 玉: 赤
maxiter = 100
result = []
for _ in 1:maxiter
x, μ, y = sample()
push!(result, y)
end
mean(result .== "赤")
0.45
maxiter = 1000000
result = []
for _ in 1:maxiter
x, μ, y = sample()
push!(result, y)
end
mean(result .== "赤")
0.399211
# 観測値（赤玉）
y_obs = "赤"
maxiter = 1_000_000
x_results = []
for _ in 1:maxiter
x, μ, y = sample()
# 生成されたyが観測と一致する場合のみ追加
y == y_obs && push!(x_results, x)
end
# 受容率（観測と一致した割合）
println("acceptance rate = $(length(x_results)/maxiter)")
# 玉が赤だった場合の袋の条件付き分布
println("p(x=A|y=赤) : approx = $(mean(x_results .== "A")) (exact=$(1/4))")
println("p(x=B|y=赤) : approx = $(mean(x_results .== "B")) (exact=$(3/4))")
acceptance rate = 0.399853 p(x=A|y=赤) : approx = 0.24998436925570147 (exact=0.25) p(x=B|y=赤) : approx = 0.7500156307442986 (exact=0.75)
function sample(N)
x = bag(rand(Bernoulli(1//2)))
μ = x=="A" ? 1//5 : 3//5
# 玉をN回復元抽出する
Y = ball.(rand(Bernoulli(μ), N))
x, μ, Y
end
sample (generic function with 2 methods)
maxiter = 10_000
Y_obs = ["赤", "赤", "白"]
x_results = []
for _ in 1:maxiter
x, μ, Y = sample(3)
# 3つの玉が完全に一致する場合のみ受容
Y == Y_obs && push!(x_results, x)
end
println("acceptance rate = $(length(x_results)/maxiter)")
println("p(x=A|y_1=赤, y_2=赤, y_3=白) : " *
"approx = $(mean(x_results .== "A")) (exact=$(2/11))")
println("p(x=B|y_1=赤, y_2=赤, y_3=白) : " *
"approx = $(mean(x_results .== "B")) (exact=$(9/11))")
acceptance rate = 0.09 p(x=A|y_1=赤, y_2=赤, y_3=白) : approx = 0.19444444444444445 (exact=0.18181818181818182) p(x=B|y_1=赤, y_2=赤, y_3=白) : approx = 0.8055555555555556 (exact=0.8181818181818182)
maxiter = 10_000
x_results = []
for _ in 1:maxiter
x, μ, Y = sample(3)
# 赤玉の個数さえ一致すれば受容するように修正
sum(Y.=="赤") == sum(Y_obs.=="赤") && push!(x_results, x)
end
println("acceptance rate = $(length(x_results)/maxiter)")
println("p(x=A|y_1=赤, y_2=赤, y_3=白) : " *
"approx = $(mean(x_results .== "A")) (exact=$(2/11))")
println("p(x=B|y_1=赤, y_2=赤, y_3=白) : " *
"approx = $(mean(x_results .== "B")) (exact=$(9/11))")
acceptance rate = 0.2638 p(x=A|y_1=赤, y_2=赤, y_3=白) : approx = 0.17930250189537528 (exact=0.18181818181818182) p(x=B|y_1=赤, y_2=赤, y_3=白) : approx = 0.8206974981046247 (exact=0.8181818181818182)
using Distributions
using PyPlot
function set_options(ax, xlabel, ylabel, title;
grid=true, gridy=false, legend=false)
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
ax.set_title(title)
if grid
if gridy
ax.grid(axis="y")
else
ax.grid()
end
end
legend && ax.legend()
end
set_options (generic function with 1 method)
# 分布を作る
d = Bernoulli(0.3)
Bernoulli{Float64}(p=0.3)
x = rand(d)
x
true
println(x)
println(Int(x))
true 1
# 10個の独立なサンプルを得る
X = rand(d, 10)
X'
1×10 adjoint(::Vector{Bool}) with eltype Bool:
0 1 1 0 0 0 0 0 1 0
# 1が生成される確率
println(pdf(d, 1))
# 0が生成する確率
println(pdf(d, 0))
# -1が生成される確率（起こり得ないのでゼロになる）
println(pdf(d, -1))
0.3 0.7 0.0
fig, ax = subplots()
ax.bar([0, 1], pdf.(d, [0,1]))
set_options(ax, "x", "probability", "probability mass function";
gridy=true)
false
println("mean = $(mean(d)), std = $(std(d))")
mean = 0.3, std = 0.458257569495584
μ = 0.3
println("std = $(sqrt(μ*(1-μ)))")
std = 0.458257569495584
X = rand(d, 10_000)
println("mean ∼ $(mean(X)), std ∼ $(std(X))")
mean ∼ 0.3019, std ∼ 0.4591050726650431
d = Binomial(20, 0.3)
Binomial{Float64}(n=20, p=0.3)
x = rand(d)
6
X = rand(d, 100)
fig, ax = subplots()
ax.hist(X)
set_options(ax, "x", "freqency", "histogram";
gridy=true)
false
println("mean (exact) = $(mean(d)), mean (approx) = $(mean(X))")
mean (exact) = 6.0, mean (approx) = 6.04
println("std (exact) = $(std(d)), std (approx) = $(std(X))")
std (exact) = 2.0493901531919194, std (approx) = 1.9741767239469057
xs = range(0, 20, length=21)
fig, ax = subplots()
ax.bar(xs, pdf.(d, xs))
set_options(ax, "x", "freqency", "probability mass function";
gridy=true)
false
# 平均パラメータのリスト
μs = [0.2, 0.4, 0.6, 0.8]
fig, axes = subplots(2, 2, sharey=true, figsize=(8,6))
for (i, ax) in enumerate(axes)
μ = μs[i]
d = Binomial(20, μ)
ax.bar(xs, pdf.(d, xs))
set_options(ax, "x", "probability", "μ=$(μ)";
gridy=true)
end
tight_layout()
M = 10
d = Multinomial(M, [0.5, 0.3, 0.2])
Multinomial{Float64, Vector{Float64}}(n=10, p=[0.5, 0.3, 0.2])
x = rand(d)
x
3-element Vector{Int64}:
6
0
4
X = rand(d, 100)
X
3×100 Matrix{Int64}:
5 5 6 3 5 3 6 4 4 4 3 4 4 … 5 3 5 5 5 3 4 4 5 6 6 4
5 3 4 4 4 2 3 5 2 2 3 3 5 5 4 3 4 3 6 3 2 3 3 3 1
0 2 0 3 1 5 1 1 4 4 4 3 1 0 3 2 1 2 1 3 4 2 1 1 5
mean(X, dims=2)
3×1 Matrix{Float64}:
4.59
3.45
1.96
mean(d)
3-element Vector{Float64}:
5.0
3.0
2.0
cov(X, dims=2)
3×3 Matrix{Float64}:
1.98172 -1.15707 -0.824646
-1.15707 1.92677 -0.769697
-0.824646 -0.769697 1.59434
cov(d)
3×3 Matrix{Float64}:
2.5 -1.5 -1.0
-1.5 2.1 -0.6
-1.0 -0.6 1.6
# 各次元の値をとる範囲
xs = 0:M
fig, ax = subplots()
cs = ax.imshow([pdf(d, [x1, x2, M - (x1 + x2)]) for x1 in xs, x2 in xs]', origin="lower")
fig.colorbar(cs)
set_options(ax, "x1", "x2", "";
grid=false)
false
μ = 2.0
d = Poisson(μ)
X = rand(d, 100)
X'
1×100 adjoint(::Vector{Int64}) with eltype Int64:
1 1 7 2 2 1 1 2 2 3 1 1 4 … 1 3 1 7 2 3 2 3 2 1 0 1
max_val = maximum(X)
fig, ax = subplots()
ax.hist(X, bins=max_val+1, range=[0, max_val])
set_options(ax, "x", "freqency", "";
gridy=true)
false
println("mean (exact) = $(mean(d)), var (exact) = $(var(d))")
mean (exact) = 2.0, var (exact) = 2.0
println("mean (approx) = $(mean(X)), var (approx) = $(var(X))")
mean (approx) = 2.05, var (approx) = 2.0277777777777777
# 表示範囲は0から25までとする
xs = 0:25
# 平均パラメータのリスト
μs = [0.5, 1.0, 4.0, 10.0]
fig, axes = subplots(2, 2, sharey=true, figsize=(8,6))
for (i, ax) in enumerate(axes)
μ = μs[i]
d = Poisson(μ)
ax.bar(xs, pdf.(d, xs))
set_options(ax, "x", "probability", "μ=$(μ)";
gridy=true)
end
tight_layout()
# 負の二項分布の作成
r = 10
μ = 0.3
d = NegativeBinomial(r, μ)
#値のサンプリング
X = rand(d, 100)
X'
1×100 adjoint(::Vector{Int64}) with eltype Int64:
13 33 17 34 11 29 20 36 18 … 14 20 17 24 22 24 13 24 18
println("mean (exact) = $(mean(d)), var (exact) = $(var(d))")
mean (exact) = 23.333333333333336, var (exact) = 77.77777777777779
println("maen(approx) = $(mean(X)), var (exact) = $(var(X))")
maen(approx) = 22.68, var (exact) = 67.81575757575757
# 表示範囲は0から60までとする
xs = 0:60
# パラメータのリスト
rs = [3, 5, 10]
μs = [0.3, 0.5, 0.7]
fig, axes = subplots(3, 3, sharey=true, figsize=(12,12))
for (i, r) in enumerate(rs)
for (j, μ) in enumerate(μs)
d = NegativeBinomial(r, μ)
axes[i, j].bar(xs, pdf.(d, xs))
# 平均と分散の計算、表示は小数3桁に丸める
m = round(mean(d), digits=3)
v = round(var(d), digits=3)
axes[i, j].text(20, 0.3, "mean=$(m), var=$(v)")
set_options(axes[i, j], "x", "probability", "r=$(r), μ=$(μ)";
gridy=true)
end
end
tight_layout()
# パラメータ
a = 0
b = 1
# 一様分布の作成
d = Uniform(a, b)
# サンプリング
X = rand(d, 1000)
X'
1×1000 adjoint(::Vector{Float64}) with eltype Float64:
0.0236739 0.252815 0.664046 0.892596 … 0.338695 0.00607204 0.631941
fig, ax = subplots()
ax.hist(X)
set_options(ax, "x", "frequency", "histogram";
gridy=true)
false
# 平均値パラメータ
μ = 0.0
# 標準偏差パラメータ
σ = 1.0
# 正規分布の作成
d = Normal(μ, σ)
# サンプリング
X = rand(d, 10000)
X'
1×10000 adjoint(::Vector{Float64}) with eltype Float64:
1.42876 -1.34033 0.559439 -1.42868 … 0.162162 -0.526015 2.3709
fig, ax = subplots()
ax.hist(X)
set_options(ax, "x", "frequency", "histogram";
gridy=true)
false
println("mean (exact) = $(mean(d)), std (exact) = $(std(d))")
println("mean (approx) = $(mean(X)), std (approx) = $(std(X))")
mean (exact) = 0.0, std (exact) = 1.0 mean (approx) = -0.009697113716219372, std (approx) = 0.9908933036261273
# 表示の範囲は−4から4までとする
xs = range(-4, 4, length=100)
# 平均パラメータのリスト
μs = [-1.0, 0.0, 1.0]
# 標準偏差パラメータリスト
σs = [0.3, 1.0, 1.5]
fig, axes = subplots(length(μs), length(σs), sharey = true, figsize=(8, 8))
for (i, μ) in enumerate(μs)
for (j, σ) in enumerate(σs)
d = Normal(μ, σ)
axes[i, j].plot(xs, pdf.(d, xs))
set_options(axes[i, j], "x", "density", "μ = $(μ), σ = $(σ)")
end
end
tight_layout()
μ = 0.0
σ = 0.2
d = Normal(μ, σ)
pdf(d, 0.1)
1.7603266338214976
xs = range(-1, 1, length=100)
fig, axes = subplots(2, 1, figsize=(4,8))
# 正規分布の確率密度関数をプロット
axes[1].plot(xs, pdf.(d, xs))
set_options(axes[1], "x", "density", "")
# 累積分布関数をプロット
axes[2].plot(xs, cdf.(d, xs))
set_options(axes[2], "x", "cumulation", "")
tight_layout()
cdf(d, 0.2) - cdf(d, 0.0)
0.34134474606854304
X = rand(d, 10000)
# ０.０から0.2にに入ったサンプルの割合を求める
mean(0.0 .< X .< 0.2)
0.3421
μ = 0.0
σ = 1.0
d = Normal(μ, σ)
Normal{Float64}(μ=0.0, σ=1.0)
X_obs = [0.1, -0.1, 0.2, 0.5]
pdf.(d, X_obs)
4-element Vector{Float64}:
0.3969525474770118
0.3969525474770118
0.3910426939754559
0.3520653267642995
prod(pdf.(d, X_obs))
0.021693249867975634
lp = sum(logpdf.(d, X_obs))
println("logpdf = $(lp)")
# 元のpdfに戻す
println("pdf = $(exp(lp))")
logpdf = -3.830754132818691 pdf = 0.021693249867975627
μ = [0.0, 0.0]
Σ = [1.0 0.0
0.0 1.0]
d = MvNormal(μ, Σ)
X = rand(d, 1000)
X
2×1000 Matrix{Float64}:
-1.57857 -0.877997 0.678171 -1.71418 … 0.293371 0.277071 2.55895
1.38254 1.35867 0.315322 0.018942 -0.87999 0.0140066 -0.510717
fig, ax = subplots(figsize=(4,4))
ax.scatter(X[1,:], X[2,:], alpha=0.5)
set_options(ax, "x1", "x2", "scatter")
false
x1s = range(-3, 3, length=100)
x2s = range(-3, 3, length=100)
fig, ax = subplots(figsize=(4,4))
cs = ax.contour(x1s, x2s, [pdf(d, [x1, x2]) for x1 in x1s, x2 in x2s]')
ax.clabel(cs, inline=true)
set_options(ax, "x1", "x2", "density")
false
μ = [0.0, 0.0]
Σ = [1.0 0.5
0.5 1.0]
d = MvNormal(μ, Σ)
FullNormal( dim: 2 μ: [0.0, 0.0] Σ: [1.0 0.5; 0.5 1.0] )
X = rand(d, 1000)
fig, ax = subplots(figsize=(4,4))
ax.scatter(X[1,:], X[2,:], alpha=0.5)
set_options(ax, "x1", "x2", "scatter")
false
x1s = range(-3, 3, length=100)
x2s = range(-3, 3, length=100)
fig, ax = subplots(figsize=(4,4))
cs = ax.contour(x1s, x2s, [pdf(d, [x1, x2]) for x1 in x1s, x2 in x2s]')
ax.clabel(cs, inline=true)
set_options(ax, "x1", "x2", "density")
false
μ = [0.0, 0.0]
Σ = [1.0 -0.5
-0.5 1.0]
d = MvNormal(μ, Σ)
FullNormal( dim: 2 μ: [0.0, 0.0] Σ: [1.0 -0.5; -0.5 1.0] )
X = rand(d, 1000)
fig, ax = subplots(figsize=(4,4))
ax.scatter(X[1,:], X[2,:], alpha=0.5)
set_options(ax, "x1", "x2", "scatter")
false
x1s = range(-3, 3, length=100)
x2s = range(-3, 3, length=100)
fig, ax = subplots(figsize=(4,4))
cs = ax.contour(x1s, x2s, [pdf(d, [x1, x2]) for x1 in x1s, x2 in x2s]')
ax.clabel(cs, inline=true)
set_options(ax, "x1", "x2", "density")
false
μ = [0.0, 0.0]
Σ = [1.5 0.25
0.25 0.5]
d = MvNormal(μ, Σ)
FullNormal( dim: 2 μ: [0.0, 0.0] Σ: [1.5 0.25; 0.25 0.5] )
function approx_integration(x_range, f)
# 幅
Δ = x_range[2] - x_range[1]
# 近似された面積と幅を返す
sum([f(x) * Δ for x in x_range]), Δ
end
approx_integration (generic function with 1 method)
# 同時分布
fig, axes = subplots(2, 2, figsize=(8,8))
cs = axes[1,1].contour(x1s, x2s, [pdf(d, [x1, x2]) for x1 in x1s, x2 in x2s]')
axes[1,1].clabel(cs, inline=true)
set_options(axes[1,1], "x1", "x2", "p(x1, x2)")
# x１の周辺分布
x_range = range(-3, 3, length=100)
p1_marginal(x1) = approx_integration(x_range, x2 -> pdf(d, [x1, x2]))[1]
axes[2,1].plot(x1s, p1_marginal.(x1s))
axes[2,1].set_xlim([minimum(x1s), maximum(x1s)])
set_options(axes[2,1], "x1", "density", "p(x1)")
# x2の周辺分布
x_range = range(-3, 3, length=100)
p2_marginal(x2) = approx_integration(x_range, x1 -> pdf(d, [x1, x2]))[1]
axes[1,2].plot(p2_marginal.(x2s), x2s)
axes[1,2].set_ylim([minimum(x2s), maximum(x2s)])
set_options(axes[1,2], "density", "x2","p(x2)")
axes[2,2].axis("off")
tight_layout()
# 同時分布
fig, axes = subplots(2, 2, figsize=(8,8))
cs = axes[1,1].contour(x1s, x2s, [pdf(d, [x1, x2]) for x1 in x1s, x2 in x2s]')
axes[1,1].clabel(cs, inline=true)
set_options(axes[1,1], "x1", "x2", "p(x1, x2)")
# x１の条件付き分布
x2 = 1.0
p1_conditional(x1) = pdf(d, [x1, x2]) / p2_marginal(x2)
axes[1,1].plot([minimum(x1s), maximum(x1s)], [x2, x2], "r--")
axes[2,1].plot(x1s, p1_conditional.(x1s), "r", label="p(x1|x2=$(x2))")
x2 = 0.0
p1_conditional(x1) = pdf(d, [x1, x2]) / p2_marginal(x2)
axes[1,1].plot([minimum(x1s), maximum(x1s)], [x2, x2], "g--")
axes[2,1].plot(x1s, p1_conditional.(x1s), "g", label="p(x1|x2=$(x2))")
axes[2,1].set_xlim([minimum(x1s), maximum(x1s)],)
set_options(axes[2,1], "x1", "density", "p(x1|x2)")
# x2の条件付き分布
x1 = 2.0
p2_conditional(x2) = pdf(d, [x1, x2]) / p1_marginal(x1)
axes[1,1].plot([x1,x1], [minimum(x2s), maximum(x2s)], "r--")
axes[1,2].plot(p2_conditional.(x2s), x2s, "r", label="p(x2|x1=$(x1))")
x1 = -2.0
p2_conditional(x2) = pdf(d, [x1, x2]) / p1_marginal(x1)
axes[1,1].plot([x1,x1], [minimum(x2s), maximum(x2s)], "g--")
axes[1,2].plot(p2_conditional.(x2s), x2s, "g", label="p(x2|x1=$(x1))")
axes[1,2].set_ylim([minimum(x2s), maximum(x2s)],)
set_options(axes[1,2], "x2", "density", "p(x2|x1)")
axes[2,2].axis("off")
tight_layout()
# ガンマ分布を作成
α = 1.5
θ = 2.5
d = Gamma(α, θ)
# サンプルの生成
X = rand(d, 100)
X'
1×100 adjoint(::Vector{Float64}) with eltype Float64:
0.813145 2.97026 7.63536 0.301381 … 4.02372 7.50194 2.93235 4.89421
fig, ax = subplots()
ax.hist(X)
set_options(ax, "x", "frequency", "histogram")
false
println("mean (exact) : $(mean(d)), $(α*θ)")
println("var (exact) : $(var(d)), $(α*θ^2)")
mean (exact) : 3.75, 3.75 var (exact) : 9.375, 9.375
xs = range(0, 10, length=100)
αs = [0.5, 1.0, 2.0]
θs = [0.5, 1.0, 1.5]
fig, axes = subplots(length(αs), length(θs), sharey=true, figsize=(8,8))
for (i, α) in enumerate(αs)
for (j, θ) in enumerate(θs)
d = Gamma(α, θ)
μ, σ = mean(d), std(d)
axes[i,j].plot(xs, pdf.(d, xs))
axes[i,j].set_ylim([0, 1.5])
set_options(axes[i,j], "x", "density", "α=$(α), θ=$(θ), \n" *
"(μ=$(round(μ, digits=3)), σ=$(round(σ, digits=3)))")
end
end
tight_layout()
# ベータ分布の作成
α = 0.5
β = 0.5
d = Beta(α, β)
# サンプルの生成
X = rand(d, 100)
X'
1×100 adjoint(::Vector{Float64}) with eltype Float64:
0.610743 0.745976 0.998187 0.525525 … 0.0209475 0.315101 0.739733
fig, ax = subplots()
ax.hist(X)
set_options(ax, "x", "frequency", "histogram")
false
xs = range(0, 1, length=100)
# パラメータのリスト
αs = [0.1, 1.0, 2.0]
βs = [0.1, 1.0, 2.0]
fig, axes = subplots(length(αs), length(βs), sharey=true, figsize=(8,8))
for (i, α) in enumerate(αs)
for (j, β) in enumerate(βs)
d = Beta(α, β)
μ, σ = mean(d), std(d)
axes[i,j].plot(xs, pdf.(d, xs))
set_options(axes[i,j], "x", "density", "α=$(α), β=$(β), \n" *
"(μ=$(round(μ, digits=3)), σ=$(round(σ, digits=3)))")
end
end
tight_layout()
# ディリクレ分布の作成
α = [0.75, 0.75, 0.75]
d = Dirichlet(α)
# サンプルの生成
X = rand(d, 1000)
X
3×1000 Matrix{Float64}:
0.280472 0.422914 0.108769 0.423537 … 0.176149 0.785769 0.0221444
0.145596 0.457917 0.435515 0.160961 0.304625 0.0926957 0.977021
0.573932 0.119168 0.455715 0.415503 0.519226 0.121535 0.000834433
fig, axes = subplots(1, 2, figsize=(8,4))
# 散布図による可視化
axes[1].scatter(X[1,:], X[2,:], alpha=0.25)
set_options(axes[1], "x1", "x2", "scatter")
# 確率密度関数
x1s = range(0, 1, length=100)
x2s = range(0, 1, length=100)
cs = axes[2].contour(x1s, x2s, [x1 + x2 > 1 ? 0.0 : pdf(d, [x1, x2, 1 - (x1 + x2)]) for x1 in x1s, x2 in x2s]')
axes[2].clabel(cs, inline=true)
set_options(axes[2], "x1", "x2", "density")
tight_layout()
mean(d)
3-element Vector{Float64}:
0.3333333333333333
0.3333333333333333
0.3333333333333333
cov(d)
3×3 Matrix{Float64}:
0.0683761 -0.034188 -0.034188
-0.034188 0.0683761 -0.034188
-0.034188 -0.034188 0.0683761
# パラメータのリスト
αs = [[0.1, 0.1, 0.1],
[0.5, 0.5, 0.5],
[1.0, 1.0, 1.0],
[2.0, 2.0, 2.0],
[5.0, 5.0, 5.0],
[0.1, 0.1, 0.5],
[0.1, 0.5, 0.1],
[0.1, 0.5, 5.0],
[1.0, 2.0, 5.0]]
xs = range(0.1, 0.99, length=100)
ys = range(0.1, 0.99, length=100)
fig, axes = subplots(3, 3, figsize=(9,9))
for (i, α) in enumerate(αs)
print(α)
d = Dirichlet(α)
X = rand(d, 1000)
axes[i].scatter(X[1,:], X[2,:], alpha=0.25)
axes[i].set_xlim([0,1])
axes[i].set_ylim([0,1])
set_options(axes[i], "x1", "x2", "α=$(α)")
end
tight_layout()
[0.1, 0.1, 0.1][0.5, 0.5, 0.5][1.0, 1.0, 1.0][2.0, 2.0, 2.0][5.0, 5.0, 5.0][0.1, 0.1, 0.5]
[0.1, 0.5, 0.1][0.1, 0.5, 5.0][1.0, 2.0, 5.0]
μ = 0.0
σ = 1.0
d = Normal(μ, σ)
X = rand(d, 100)
Y = exp.(X)
Y'
1×100 adjoint(::Vector{Float64}) with eltype Float64:
0.856311 0.832838 3.73996 4.16099 … 2.29158 0.986543 1.66065 4.80166
fig, axes = subplots(1, 2, figsize=(12,4))
axes[1].hist(X)
set_options(axes[1], "x", "freqency", "Normal distribution";
gridy=true)
axes[2].hist(Y)
set_options(axes[2], "y", "frequency", "Log-Normal distribution";
gridy=true)
tight_layout()
d = LogNormal(μ, σ)
X = rand(d, 100)
X'
1×100 adjoint(::Vector{Float64}) with eltype Float64:
0.558211 2.9403 0.564984 0.803061 … 1.96972 1.23392 0.1873 2.36384
fig, ax = subplots()
ax.hist(X)
set_options(ax, "x", "freqency", "histogram";
gridy=true)
false
xs = range(0, 5, length=100)
μs = [-1, 0.0, 1.0]
σs = [0.2, 1.0, 1.5]
fig, axes = subplots(length(μs), length(σs), sharey=true, figsize=(12,12))
for (i, μ) in enumerate(μs)
for (j, σ) in enumerate(σs)
d = LogNormal(μ, σ)
axes[i,j].plot(xs, pdf.(d, xs))
m = round(mean(d), digits=2)
v = round(var(d), digits=2)
axes[i,j].text(2, 5, "mean = $(m), var = $(v)")
set_options(axes[i,j], "x", "density", "μ = $(μ), σ = $(σ)")
end
end
tight_layout()
# パラメータを生成するための分布
μ1 = [-1.0, 1.0]
Σ1 = [0.2 0.0
0.0 0.2]
μ2 = [1.0, -1.0]
Σ2 = [0.4 0.0
0.0 0.4]
p = 0.3
# サンプリングする回数
N = 100
# 生成データを保存する配列
X = Array{Float64}(undef, 2, N)
# 選択された分布を示す配列（潜在変数）
S = Array{Bool}(undef, N)
for i in 1:N
# 潜在変数のサンプル
S[i] = rand(Bernoulli(p))
# 潜在変数の値に応じて多変量正規分布のパラメータを切り替える
(μ, Σ) = S[i] == 1 ? (μ1, Σ1) : (μ2, Σ2)
# データをサンプリングする
X[:, i] = rand(MvNormal(μ, Σ))
end
# 生成されたデータを散布図として可視化
fig, ax = subplots(figsize=(4,4))
ax.scatter(X[1,:], X[2,:], alpha=0.5)
set_options(ax, "x1", "x2", "scatter")
false
# 混合多変量正規分布の確率密度関数
pdfgmm(x) = p*pdf(MvNormal(μ1, Σ1), x) + (1-p)*pdf(MvNormal(μ2, Σ2), x)
xs1 = range(-3, 3, length=100)
xs2 = range(-3, 3, length=100)
fig, ax = subplots(figsize=(4,4))
cs = ax.contour(xs1, xs2, [pdfgmm([x1, x2]) for x1 in xs1, x2 in xs2]')
ax.clabel(cs, inline=true)
set_options(ax, "x1", "x2", "density")
false
# パラメータを生成するための分布
μ = [0.0, 0.0]
Σ = [0.1 0.0
0.0 0.1]
# 出力値に付加するノイズ
σ = 1.0
# 入力値はあらかじめ与えておく
X = [-10, -5, 0, 5, 10]
# サンプリングする回数
num_samples = 3
# パラメータのサンプル
W = rand(MvNormal(μ, Σ), num_sample)
# 出力値を保存するためのリスト
Ys = []
fig, axes = subplots(1, 2, figsize=(8,3))
xs = range(-12, 12, length=100)
for n in 1:num_samples
w1, w2 = W[:,n]
# パラメータをプロット
axes[1].scatter(w1, w2)
# 生成された関数をプロット
linear_function(x) = w1*x + w2
axes[2].plot(xs, linear_function.(xs))
# 関数からの出力値も生成する
Y = rand.(Normal.(linear_function.(X), σ))
push!(Ys, Y)
end
set_options(axes[1], "w1", "w2", "parameters")
set_options(axes[2], "x", "y", "functions")
tight_layout()
UndefVarError: num_sample not defined Stacktrace: [1] top-level scope @ In[207]:16 [2] eval @ ./boot.jl:360 [inlined] [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String) @ Base ./loading.jl:1116
fig, axes = subplots(1, num_sample, sharey=true, figsize=(12,3))
for n in 1:num_sample
w1, w2 = W[:,n]
Y = Ys[n]
linear_function(x) =w1*x + w2
axes[n].plot(xs, linear_function.(xs))
axes[n].scatter(X, Y)
set_options(axes[n], "x", "y", "w1 = $(round(w1, digits=3)), w2 = $(round(w2, digits=3))")
end
tight_layout()
UndefVarError: num_sample not defined Stacktrace: [1] top-level scope @ In[208]:1 [2] eval @ ./boot.jl:360 [inlined] [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String) @ Base ./loading.jl:1116
sig(x) = 1 / (1 + exp(-x))
sig (generic function with 1 method)
# パラメータを生成するための分布
μ = [0.0, 0.0]
Σ = [0.01 0.0
0.0 0.01]
# 入力値はあらかじめ与えておく
X = [-10, -5, 0, 5, 10]
# サンプリングする回数
num_samples = 3
# パラメータのサンプル
W = rand(MvNormal(μ, Σ), num_sample)
# 出力値を保存するためのリスト
Ys = []
fig, axes = subplots(1, 2, figsize=(8,3))
xs = range(-12, 12, length=100)
for n in 1:num_samples
w1, w2 = W[:,n]
# パラメータをプロット
axes[1].scatter(w1, w2)
# 生成された関数をプロット
nonlinear_function(x) = sig(w1*x + w2)
axes[2].plot(xs, nonlinear_function.(xs))
# 関数からの出力値も生成する
Y = rand.(Bernoulli.(nonlinear_function.(X)))
push!(Ys, Y)
end
set_options(axes[1], "w1", "w2", "parameters")
set_options(axes[2], "x", "y", "functions")
tight_layout()
UndefVarError: num_sample not defined Stacktrace: [1] top-level scope @ In[210]:13 [2] eval @ ./boot.jl:360 [inlined] [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String) @ Base ./loading.jl:1116
fig, axes = subplots(1, num_sample, sharey=true, figsize=(12,3))
xs = range(-12, 12, length=100)
for n in 1:num_samples
w1, w2 = W[:,n]
Y = Ys[n]
nonlinear_function(x) = sig(w1*x + w2)
axes[n].plot(xs, nonlinear_function.(xs))
axes[n].scatter(X, Y)
set_options(axes[n], "x", "y", "w1=$(round(w1, digits=3)), w2=$(round(w2, digits=3))")
end
tight_layout()
UndefVarError: num_sample not defined Stacktrace: [1] top-level scope @ In[211]:1 [2] eval @ ./boot.jl:360 [inlined] [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String) @ Base ./loading.jl:1116
# パラメータを生成するための分布
μ = [0.0, 0.0]
Σ = [0.01 0.0
0.0 0.01]
# 入力値はあらかじめ与えておく
X = [-10, -5, 0, 5, 10]
# サンプリングする回数
num_samples = 3
# パラメータのサンプル
W = rand(MvNormal(μ, Σ), num_sample)
# 出力値を保存するためのリスト
Ys = []
fig, axes = subplots(1, 2, figsize=(8,3))
xs = range(-12, 12, length=100)
for n in 1:num_samples
w1, w2 = W[:,n]
# パラメータをプロット
axes[1].scatter(w1, w2)
# 生成された関数をプロット
nonlinear_function(x) = exp(w1*x + w2)
axes[2].plot(xs, nonlinear_function.(xs))
# 関数からの出力値も生成する
Y = rand.(Poisson.(nonlinear_function.(X)))
push!(Ys, Y)
end
set_options(axes[1], "w1", "w2", "parameters")
set_options(axes[2], "x", "y", "functions")
tight_layout()
UndefVarError: num_sample not defined Stacktrace: [1] top-level scope @ In[212]:13 [2] eval @ ./boot.jl:360 [inlined] [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String) @ Base ./loading.jl:1116
fig, axes = subplots(1, num_sample, sharey=true, figsize=(12,3))
xs = range(-12, 12, length=100)
for n in 1:num_samples
w1, w2 = W[:,n]
Y = Ys[n]
nonlinear_function(x) = exp(w1*x + w2)
axes[n].plot(xs, nonlinear_function.(xs))
axes[n].scatter(X, Y)
set_options(axes[n], "x", "y", "w1=$(round(w1, digits=3)), w2=$(round(w2, digits=3))")
end
tight_layout()
UndefVarError: num_sample not defined Stacktrace: [1] top-level scope @ In[213]:1 [2] eval @ ./boot.jl:360 [inlined] [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String) @ Base ./loading.jl:1116
# パッケージの読み込み
using Distributions, PyPlot, LinearAlgebra
WARNING: using LinearAlgebra.I in module Main conflicts with an existing identifier.
function generate(N)
μ = rand(Uniform(0, 1))
X = rand(Bernoulli(μ), N)
μ, X
end
# 5回コイン投げを行う
generate(5)
(0.14466576460977287, Bool[0, 0, 0, 1, 0])
# 1を表、０を裏とする
side(x) = x == 1 ? "表" : "裏"
for i in 1:10
μ, X = generate(5)
println("コイン$(i)、表が出る確率 μ = $(μ)、出目 X = $(side.(X))")
end
コイン1、表が出る確率 μ = 0.2115239161862532、出目 X = ["表", "表", "表", "裏", "裏"] コイン2、表が出る確率 μ = 0.19157334256084235、出目 X = ["表", "表", "裏", "表", "表"] コイン3、表が出る確率 μ = 0.6340337570025956、出目 X = ["表", "表", "表", "裏", "表"] コイン4、表が出る確率 μ = 0.4673419080511929、出目 X = ["裏", "表", "表", "表", "裏"] コイン5、表が出る確率 μ = 0.7624851169457685、出目 X = ["表", "裏", "表", "裏", "裏"] コイン6、表が出る確率 μ = 0.7451980753784295、出目 X = ["表", "裏", "表", "裏", "表"] コイン7、表が出る確率 μ = 0.26375329479765175、出目 X = ["表", "裏", "表", "裏", "表"] コイン8、表が出る確率 μ = 0.8119528312066475、出目 X = ["裏", "表", "表", "表", "表"] コイン9、表が出る確率 μ = 0.6140591793610704、出目 X = ["表", "裏", "裏", "表", "表"] コイン10、表が出る確率 μ = 0.33996595666852136、出目 X = ["裏", "表", "裏", "裏", "裏"]
μs = range(0, 1, length=100)
fig, ax = subplots()
ax.plot(μs, pdf.(Uniform(0, 1), μs))
set_options(ax, "μ", "density", "Uniform distribution")
ax.set_xlim([0, 1])
ax.set_ylim([0, 1.1])
(0.0, 1.1)
# 「裏、裏、裏、表、表」をデータとして取得
X_obs1 = [0,0,0,1,1]
5-element Vector{Int64}:
0
0
0
1
1
maxiter = 1_000_000
μ_posteriori1 = []
for i in 1:maxiter
# パラメータおよびデータの生成
μ, X = generate(length(X_obs1))
# X内の1の合計が観測と一致していれば、この時のパラメータを受容
sum(X) == sum(X_obs1) && push!(μ_posteriori1, μ)
end
# 受容率の計算
acceptance_rate = length(μ_posteriori1) / maxiter
println("acceptance rate = $(acceptance_rate)")
μ_posteriori'
acceptance rate = 0.166686
UndefVarError: μ_posteriori not defined Stacktrace: [1] top-level scope @ In[219]:15 [2] eval @ ./boot.jl:360 [inlined] [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String) @ Base ./loading.jl:1116
fig, ax = subplots()
ax.hist(μ_posteriori1)
ax.set_xlim([0,1])
set_options(ax, "μ", "frequency", "histogram";
gridy=true)
false
X_obs2 = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1]
20-element Vector{Int64}:
0
0
0
0
0
0
0
0
0
0
0
0
1
1
1
1
1
1
1
1
maxiter = 1_000_000
μ_posteriori2 = []
for i in 1:maxiter
# パラメータおよびデータの生成
μ, X = generate(length(X_obs2))
# X内の1の合計が観測と一致していれば、この時のパラメータを受容
sum(X) == sum(X_obs2) && push!(μ_posteriori2, μ)
end
# 受容率の計算
acceptance_rate = length(μ_posteriori2) / maxiter
println("acceptance rate = $(acceptance_rate)")
μ_posteriori'
acceptance rate = 0.047841
UndefVarError: μ_posteriori not defined Stacktrace: [1] top-level scope @ In[222]:15 [2] eval @ ./boot.jl:360 [inlined] [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String) @ Base ./loading.jl:1116
fig, ax = subplots()
ax.hist(μ_posteriori2)
ax.set_xlim([0,1])
set_options(ax, "μ", "frequency", "histogram";
gridy=true)
false
pred1 = mean(rand.(Bernoulli.(μ_posteriori1)))
pred2 = mean(rand.(Bernoulli.(μ_posteriori2)))
println("$(pred1), $(pred2)")
0.4280563454639262, 0.40691039066909135
function generate2(N)
# μの事前分布を修正
μ = rand(Uniform(0, 0.5))
X = rand(Bernoulli(μ), N)
μ, X
end
generate2(5)
(0.33987893554338955, Bool[0, 0, 0, 0, 1])
for i in 1:10
μ, X = generate2(5)
println("コイン $(i)、表が出る確率 μ = $(μ)、出目 X = $(side.(X))")
end
コイン 1、表が出る確率 μ = 0.371291143488243、出目 X = ["表", "裏", "裏", "表", "裏"] コイン 2、表が出る確率 μ = 0.11257772410611666、出目 X = ["裏", "裏", "裏", "表", "表"] コイン 3、表が出る確率 μ = 0.4383712405951564、出目 X = ["表", "裏", "表", "表", "裏"] コイン 4、表が出る確率 μ = 0.23543633208856873、出目 X = ["裏", "裏", "裏", "裏", "裏"] コイン 5、表が出る確率 μ = 0.0069485910271677165、出目 X = ["裏", "裏", "裏", "裏", "裏"] コイン 6、表が出る確率 μ = 0.05611836072281462、出目 X = ["裏", "裏", "裏", "裏", "裏"] コイン 7、表が出る確率 μ = 0.4629175852579448、出目 X = ["表", "表", "裏", "表", "裏"] コイン 8、表が出る確率 μ = 0.35259901883072686、出目 X = ["裏", "表", "表", "表", "裏"] コイン 9、表が出る確率 μ = 0.03320618117563745、出目 X = ["裏", "裏", "裏", "裏", "裏"] コイン 10、表が出る確率 μ = 0.15212636716313133、出目 X = ["裏", "裏", "裏", "裏", "裏"]
maxiter = 1_000_000
μ_posteriori3 = []
for i in 1:maxiter
# パラメータおよびデータの生成
μ, X = generate2(length(X_obs1))
# X内の1の合計が観測と一致していれば、この時のパラメータを受容
sum(X) == sum(X_obs1) && push!(μ_posteriori3, μ)
end
# 受容率の計算
acceptance_rate = length(μ_posteriori3) / maxiter
println("acceptance rate = $(acceptance_rate)")
μ_posteriori'
acceptance rate = 0.219677
UndefVarError: μ_posteriori not defined Stacktrace: [1] top-level scope @ In[227]:15 [2] eval @ ./boot.jl:360 [inlined] [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String) @ Base ./loading.jl:1116
fig, ax = subplots()
ax.hist(μ_posteriori3)
ax.set_xlim([0,1])
set_options(ax, "μ", "frequency", "histogram";
gridy=true)
false
# 同時分布p(X, μ)の確率密度関数
p_joint(X, μ) = prod(pdf.(Bernoulli(μ), X)) * pdf(Uniform(0,1), μ)
# 数値積分
function approx_integration(μ_range, p)
Δ = μ_range[2] - μ_range[1]
X -> sum([p(X, μ) * Δ for μ in μ_range]), Δ
end
# μの積分範囲
μ_range = range(0, 1, length=100)
# 数値積分の実行
p_marginal, Δ = approx_integration(μ_range, p_joint)
# データ（２種類）
X_obs1 = [0,0,0,1,1]
X_obs2 = [0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1]
# それぞれの周辺尤度の近似計算
println("$(p_marginal(X_obs1)), $(p_marginal(X_obs2))")
0.016666666493163274, 3.7801895387034807e-7
# パラメータの可視化範囲
μs = range(0, 1, length=100)
fig, axes = subplots(1, 2, sharey=true, figsize=(10,4))
for (i, X_obs) in enumerate([X_obs1, X_obs2])
posterior(μ) = p_joint(X_obs, μ) / p_marginal(X_obs)
axes[i].plot(μs, posterior.(μs))
set_options(axes[i], "μ", "density", "approximate posterior (X_obs$(i))")
end
# 積分の中身の式
posterior1(μ) = p_joint(X_obs1, μ) / p_marginal(X_obs1)
posterior2(μ) = p_joint(X_obs2, μ) / p_marginal(X_obs2)
p_inner1(x, μ) = pdf.(Bernoulli(μ), x) * posterior1(μ)
p_inner2(x, μ) = pdf.(Bernoulli(μ), x) * posterior2(μ)
# パラメータμに関する積分
μ_range = range(0, 1, length=100)
pred1, Δ1 = approx_integration(μ_range, p_inner1)
pred2, Δ2 = approx_integration(μ_range, p_inner2)
# 1が出る予測確率
println("$(pred1(1)), $(pred2(1))")
0.4285714434416308, 0.40909090909090784
fig, axes = subplots(1, 2, sharey=true, figsize=(10,4))
μs = range(0, 1, length=100)
for (i, X_obs) in enumerate([X_obs1, X_obs2])
# 厳密な事後分布はベータ分布
α = 1.0 + sum(X_obs)
β = 1.0 + length(X_obs) - sum(X_obs)
d = Beta(α, β)
# 事後分布を可視化
axes[i].plot(μs, pdf.(d, μs))
set_options(axes[i], "μ", "density", "exact posterior (X_obs$(i))")
end
function prediction(X_obs)
α = 1.0 + sum(X_obs)
β = 1.0 + length(X_obs) - sum(X_obs)
α / (α + β)
end
println("$(prediction(X_obs1)), $(prediction(X_obs2))")
0.42857142857142855, 0.4090909090909091
function generate_linear(X, σ, μ1, μ2, σ1, σ2)
w1 = rand(Normal(μ1, σ1))
w2 = rand(Normal(μ2, σ2))
linear_function(x) = w1*x + w2
Y = rand.(Normal.(linear_function.(X), σ))
Y, linear_function, w1, w2
end
generate_linear (generic function with 1 method)
# あらかじめ与えられるパラメータと入力集合X
σ = 1.0
μ1 = 0.0
μ2 = 0.0
σ1 = 10.0
σ2 = 10.0
X = [-1.0, -0.5, 0, 0.5, 1.0]
# 可視化する範囲
xs = range(-2, 2, length=100)
fig, axes = subplots(2, 3, sharey=true, figsize=(12,6))
for ax in axes
# 関数f、出力Yの生成
Y, linear_function, w1, w2 = generate_linear(X, σ, μ1, μ2, σ1, σ2)
# 生成された直線とYのプロット
ax.plot(xs, linear_function.(xs), label="simulated function")
ax.scatter(X, Y, label="simulated data")
set_options(ax, "x", "y", "data (N = $(length(X)))", legend=true)
end
tight_layout()
# 平均パラメータのリスト
μ1_list = [-20.0, 0.0, 20.0]
μ2_list = [-20.0, 0.0, 20.0]
# 標準偏差パラメータは固定
σ1 = 10.0
σ2 = 10.0
fig, axes = subplots(length(μ1_list), length(μ2_list), sharey=true, figsize=(12, 12))
for (i, μ1) in enumerate(μ1_list)
for(j, μ2) in enumerate(μ2_list)
# 関数を複数サンプル
fs = [generate_linear(X, σ, μ1, μ2, σ1, σ2)[2] for _ in 1:100]
# 生成された直線群をプロット
for linear_function in fs
axes[i,j].plot(xs, linear_function.(xs), "g", alpha=0.1)
end
axes[i,j].set_xlim(extrema(xs))
set_options(axes[i,j], "x", "y", "μ1=$(μ1), μ2=$(μ2)")
end
end
tight_layout()
# 標準偏差パラメータのリスト
σ1_list = [1.0, 10.0, 20.0]
σ2_list = [1.0, 10.0, 20.0]
# 平均パラメータは固定
μ1 = 0.0
μ2 = 0.0
fig, axes = subplots(length(σ1_list), length(σ2_list), sharey=true, figsize=(12, 12))
for (i, σ1) in enumerate(σ1_list)
for(j, σ2) in enumerate(σ2_list)
# 関数を複数サンプル
fs = [generate_linear(X, σ, μ1, μ2, σ1, σ2)[2] for _ in 1:100]
# 生成された直線群をプロット
for linear_function in fs
axes[i,j].plot(xs, linear_function.(xs), "g", alpha=0.1)
end
axes[i,j].set_xlim(extrema(xs))
set_options(axes[i,j], "x", "y", "σ1=$(σ1), σ2=$(σ2)")
end
end
tight_layout()
# 入力データ
X_obs = [-2, 1, 5]
# 出力データ
Y_obs = [-2.2, -1.0, 1.5]
# 散布図で可視化
fig, ax = subplots()
ax.scatter(X_obs, Y_obs)
set_options(ax, "x", "y", "data = (N = $(length(X_obs)))")
false
# 同時分布
p_joint(X, Y, w1, w2) = prod(pdf.(Normal.(w1.*X .+ w2, σ), Y)) * pdf(Normal(μ1, σ1), w1) * pdf(Normal(μ2, σ2), w2)
# 数値積分
function approx_integration_2D(w_range, p)
Δ = w_range[2] - w_range[1]
(X, Y) -> sum([p(X, Y, w1, w2) * Δ^2 for w1 in w_range, w2 in w_range])
end
# wの積分範囲
w_range = range(-3, 3, length=100)
# 数値積分の実行
p_marginal = approx_integration_2D(w_range, p_joint)
p_marginal(X_obs, Y_obs)
6.924264340150274e-5
# 事後分布の計算
w_posterior = [p_joint(X_obs, Y_obs, w1, w2) for w1 in w_range, w2 in w_range] ./ p_marginal(X_obs, Y_obs)
fig, axes = subplots(1, 2, figsize=(8,4))
# 等高線図
cs = axes[1].contour(w_range, w_range, w_posterior', cmap="jet")
axes[1].clabel(cs, inline=true)
set_options(axes[1], "w1", "w2", "posterior density (contour)")
# カラーメッシュ
xgrid = repeat(w_range', length(w_range), 1)
ygrid = repeat(w_range, 1, length(w_range))
axes[2].pcolormesh(xgrid, ygrid, w_posterior', cmap="jet", shading="auto")
set_options(axes[2], "w1", "w2", "posterior density (colored)")
tight_layout()
function approx_predictive(w_posterior, w_range, p)
Δ = w_range[2] - w_range[1]
(x, y) -> sum([p(x, y, w1, w2) * w_posterior[i,j] * Δ^2 for (i, w1) in enumerate(w_range), (j, w2) in enumerate(w_range)])
end
p_likelihood(xp, yp, w1, w2) = pdf(Normal(w1*xp + w2, σ), yp)
p_predictive = approx_predictive(w_posterior, w_range, p_likelihood)
#80 (generic function with 1 method)
xp = 4.0
fig, ax = subplots()
ys = range(-5, 5, length=100)
ax.plot(ys, p_predictive.(xp, ys))
set_options(ax, "y_p", "density", "predictive distiribution p(y_p|x_p=$(xp), X=X_obs, Y=Y_obs)")
false
# 描画範囲
xs = range(-10, 10, length=100)
ys = range(-5, 5, length=100)
# 密度の計算
density_y = p_predictive.(xs, ys')
fig, axes = subplots(1, 2, sharey=true, figsize=(8,4))
# 等高線図
cs = axes[1].contour(xs, ys, density_y', cmap="jet")
axes[1].clabel(cs, inline=true)
axes[1].scatter(X_obs, Y_obs)
set_options(axes[1], "x", "y", "predictive distribution (contour)")
# カラーメッシュ
xgrid = repeat(xs', length(ys), 1)
ygrid = repeat(ys, 1, length(xs))
axes[2].pcolormesh(xgrid, ygrid, density_y', cmap="jet", shading="auto")
axes[2].plot(X_obs, Y_obs, "ko", label="data")
set_options(axes[2], "x", "y", "predictive distribution (colored)")
tight_layout()
# 切片用の擬似データ
X_extended = hcat(X_obs, ones(size(X_obs)))'
# 事前分布のパラメータ
Σ = [σ1^2 0
0 σ2^2]
μ = [μ1, μ2]
# 事後分布のパラメータ
Σ_hat = inv(σ^(-2) * X_extended*X_extended' + inv(Σ))
μ_hat = Σ_hat*(σ^(-2) * X_extended*Y_obs + inv(Σ)*μ)
# 予測分布のパラメータ
μp(xp) = (μ_hat' * xp)[1]
σp(xp) = sqrt(σ^2 + (xp' * Σ_hat * xp)[1])
# 予測分布の確率密度関数
p_predictive_exact(xp, yp) = pdf(Normal(μp(xp), σp(xp)), yp)
# 描画範囲
xs = range(-10, 10, length=100)
ys = range(-5, 5, length=100)
density_y = [p_predictive_exact([x, 1.0], y) for x in xs, y in ys]
fig, axes = subplots(1, 2, sharey=true, figsize=(8,4))
# 等高線図
cs = axes[1].contour(xs, ys, density_y', cmap="jet")
axes[1].clabel(cs, inline=true)
axes[1].scatter(X_obs, Y_obs)
set_options(axes[1], "x", "y", "predictive distiribution (contour)")
# カラーメッシュ
xgrid = repeat(xs', length(ys), 1)
ygrid = repeat(ys, 1, length(xs))
axes[2].pcolormesh(xgrid, ygrid, density_y', cmap="jet", shading="auto")
axes[2].plot(X_obs, Y_obs, "ko", label="data")
set_options(axes[2], "x", "y", "predictive distribution (colored)")
tight_layout()
# シグモイド関数を定義
sig(x) = 1/(1+exp(-x))
# パラメータ、関数、出力集合Yを生成
function generate_logistic(X, μ1, μ2, σ1, σ2)
w1 = rand(Normal(μ1, σ1))
w2 = rand(Normal(μ2, σ2))
nonlinear_function(x) = sig(w1*x + w2)
Y = rand.(Bernoulli.(nonlinear_function.(X)))
Y, nonlinear_function, w1, w2
end
generate_logistic (generic function with 1 method)
# あらかじめ与えられたパラメータと入力X
μ1 = 0
μ2 = 0
σ1 = 10.0
σ2 = 10.0
X = [-1.0, -0.5, 0, 0.5, 1.0]
# 可視化する範囲
xs = range(-2, 2, length=100)
fig, axes = subplots(2, 3, sharey=true, figsize=(12,6))
for ax in axes
# 関数f、出力Yの生成
Y, nonlinear_function, w1, w2 = generate_logistic(X, μ1, μ2, σ1, σ2)
# 生成された直線とYのプロット
ax.plot(xs, nonlinear_function.(xs), label="simulated function")
ax.scatter(X, Y, label="simulated data")
set_options(ax, "x", "y", "data (N = $(length(X)))", legend=true)
end
tight_layout()
# 入力データセット
X_obs = [-2, 1, 2]
# 出力データセット
Y_obs = Bool.([0, 1, 1])
# 散布図で可視化
fig, ax = subplots()
ax.scatter(X_obs, Y_obs)
ax.set_xlim([-3, 3])
set_options(ax, "x", "y", "data (N = $(length(X_obs)))")
false
# 最大サンプリング数
maxiter = 10_000
# パラメータ保存用
param_posterior = Vector{Tuple{Float64, Float64}}()
for i in 1:maxiter
# 関数f、出力Yの生成
Y, nonlinear_function, w1, w2 = generate_logistic(X_obs, μ1, μ2, σ1, σ2)
# 観測データと一致していれば、そのときのパラメータwを保存
Y == Y_obs && push!(param_posterior, (w1, w2))
end
# サンプル受容率
acceptance_rate = length(param_posterior) / maxiter
println("acceptance_rate = $(acceptance_rate)")
acceptance_rate = 0.2979
# パラメータ抽出用の関数
unzip(a) = map(x -> getfield.(a, x), fieldnames(eltype(a)))
# 事前分布からのサンプル（10,000個）
param_prior = [generate_logistic(X, μ1, μ2, σ1, σ2)[3:4] for i in 1:10_000]
w1_prior, w2_prior = unzip(param_prior)
# 事後分布からのサンプル
w1_posterior, w2_posterior = unzip(param_posterior)
fig, axes = subplots(1, 2, sharex=true, sharey=true, figsize=(8,4))
# 事前分布
axes[1].scatter(w1_prior, w2_prior, alpha=0.01)
set_options(axes[1], "w1", "w2", "sample from prior")
# 事後分布
axes[2].scatter(w1_posterior, w2_posterior, alpha=0.01)
set_options(axes[2], "w1", "w2", "sample from posterior")
tight_layout()
# 関数を可視化する範囲
xs = range(-3, 3, length=100)
# サンプリングで得られたパラメータ全体のプロット
fig, ax = subplots()
ax.scatter(w1_posterior, w2_posterior, alpha=0.1)
set_options(ax, "w1", "w2", "sample from posterior")
fig, axes = subplots(2, 3, figsize=(12,8))
for i in eachindex(axes)
# 関数を可視化するためのwを1つ適当に選択
j = round(Int, length(param_posterior)*rand()) + 1
w1, w2 = param_posterior[j]
# 選択されたw
ax.scatter(w1, w2, color="r")
ax.text(w1, w2, i)
# 対応する関数のプロット
nonlinear_function(x) = sig(w1*x + w2)
axes[i].plot(xs, nonlinear_function.(xs), "r")
# 観測データのプロット
axes[i].scatter(X_obs, Y_obs)
axes[i].set_ylim([-0.1, 1.1])
set_options(axes[i], "x", "y (prob.)", "($(i)) w1=$(round(w1, digits=3)), w2=$(round(w2, digits=3))")
end
tight_layout()
fig, axes = subplots(1, 2, sharey=true, figsize=(12,4))
fs = []
for (i, param) in enumerate(param_posterior)
# 1サンプル分のパラメータ
w1, w2 = param
# 1サンプル分の予測関数
nonlinear_function(x) = sig(w1*x + w2)
axes[1].plot(xs, nonlinear_function.(xs), "g", alpha=0.01)
push!(fs, nonlinear_function.(xs))
end
axes[1].scatter(X_obs, Y_obs)
set_options(axes[1], "x", "y (prob.)", "function samples")
# 予測確率
axes[2].plot(xs, mean(fs), label="prediction")
axes[2].scatter(X_obs, Y_obs, label="data")
set_options(axes[2], "x", "y (prob.)", "prediction", legend=true)
tight_layout()
# 予測対象の点候補
x_list = [-1.0, 0.0, 1.0]
fig_num = length(x_list)
fig, axes = subplots(fig_num, 2, sharey=true, figsize=(12, 4*fig_num))
for (j, x) in enumerate(x_list)
# パラメータごとに関数を可視化
for (i, param) in enumerate(param_posterior)
w1, w2 = param
nonlinear_function(x) = sig(w1*x + w2)
axes[j].plot(xs, nonlinear_function.(xs), "g", alpha=0.01)
end
# 観測データ
axes[j].scatter(X_obs, Y_obs, label="data")
# 候補点のx座標
axes[j].plot([x, x], [0, 1], "r--", label="x_p=$(x)")
axes[j].set_xlim(extrema(xs))
set_options(axes[j], "x", "y (prob.)", "function sample")
axes[j].legend(loc="lower right")
# 点xにおける関数値（確率値）をヒストグラムとして可視化
probs = [sig(param[1]*x+param[2]) for param in param_posterior]
axes[j+fig_num].hist(probs, orientation="horizontal")
set_options(axes[j+fig_num], "y", "freqency", "function samples";
gridy=true)
axes[j+fig_num].grid(axis="x")
end
# 同時分布p(Y, w|X)
p_joint(X, Y, w1, w2) = prod(pdf.(Bernoulli.(sig.(w1.*X .+ w2)), Y)) * pdf(Normal(μ1, σ1), w1) * pdf(Normal(μ2, σ2), w2)
# μの積分範囲
w_range = range(-30, 30, length=100)
# 数値積分の実行
p_marginal = approx_integration_2D(w_range, p_joint)
p_marginal(X_obs, Y_obs)
0.2966829565005522
w_posterior = [p_joint(X_obs, Y_obs, w1, w2) for w1 in w_range, w2 in w_range] ./ p_marginal(X_obs, Y_obs)
# 事後分布の描画
fig, axes = subplots(1, 2, figsize=(8,4))
cs = axes[1].contour(w_range, w_range, w_posterior' .+ eps(), cmap="jet")
axes[1].clabel(cs, inline=true)
set_options(axes[1], "w1", "w2", "posterior density (contour)")
xgrid = repeat(w_range', length(w_range), 1)
ygrid = repeat(w_range, 1, length(w_range))
axes[2].pcolormesh(xgrid, ygrid, w_posterior', cmap="jet", shading="auto")
set_options(axes[2], "w1", "w2", "posteior density (colored)")
tight_layout()
p_likelihood(xp, yp, w1, w2) = pdf(Bernoulli(sig(w1*xp + w2)), yp)
p_predictive = approx_predictive(w_posterior, w_range, p_likelihood)
#80 (generic function with 1 method)
# 関数を可視化する範囲
xs = range(-3, 3, length=50)
fig, ax = subplots()
ax.scatter(X_obs, Y_obs, label="data")
ax.plot(xs, p_predictive.(xs, 1), label="probability")
ax.set_xlim([-3, 3])
set_options(ax, "x", "y", "prediction", legend=true)
PyObject <matplotlib.legend.Legend object at 0x7fe05b723ac0>
using Distributions, PyPlot, LinearAlgebra
# n次元単位行列
# eye(n) = Diagonal{Float64}(I, n)
eye(n) = Diagonal{Float64}(ones(n))
# パラメータ抽出用の関数
unzip(a) = map(a -> getfield.(a,x), fieldnames(eltype(a)))
unzip (generic function with 1 method)
# 入力データセット
X_obs = [-2, 1, 5]
# 出力データセット
Y_obs = [-2.2, -1.0, 1.5]
# 散布図で可視化
fig, ax = subplots()
ax.scatter(X_obs, Y_obs)
set_options(ax, "x", "y", "data (N = $(length(X_obs)))")
false
# 切片はゼロで固定
w2 = 0
# yに付加されているノイズの標準偏差
σ = 1.0
# 事前分布の平均値と標準偏差
μ1 = 0.0
σ1 = 10.0
10.0
# 非正規化対数事後分布
ulp(w1) = sum(logpdf.(Normal.(w1*X_obs .+ w2, σ), Y_obs)) + logpdf(Normal(μ1, σ1), w1)
# 分布を表示する範囲
w1s = range(-5, 5, length=100)
fig, axes = subplots(1, 2, figsize=(8,4))
# 非正規化対数事後分布の可視化
axes[1].plot(w1s, ulp.(w1s))
set_options(axes[1], "w1", "log density (unnormalized)", "unnormarized log posterior distribution")
# 非正規化事後分布の可視化
axes[2].plot(w1s, exp.(ulp.(w1s)))
set_options(axes[2], "w1", "density (unnormalized)", "unnormarized posterior distribution")
tight_layout()
# 勾配法（1dim）
function gradient_method_1dim(funct, x_init, η, maxiter)
ffunct(x) = ForwardDiff.derivative(funct, x)
x_seq = Array{typeof(x_init), 1}(undef, maxiter)
x_seq[1] = x_init
for i in 2:maxiter
x_seq[i] = x_seq[i-1] + η*ffunct(x_seq[i-1])
end
x_seq
end
# 最適化パラメータ
w1_init = 0.0
maxiter = 100
η = 0.01
# 最適化の実施
w1_seq = gradient_method_1dim(ulp, w1_init, η, maxiter)
# 勾配法の過程を可視化
fig, ax = subplots(figsize=(8,4))
ax.plot(w1_seq)
set_options(ax, "iteration", "w1", "w1 sequence")
false
# 最適化パラメータ
w1_init = 0.0
maxiter = 100
η = 0.01
# 最適化のラッパー関数の定義
function inference_wrapper_gd_1dim(log_joint, params, w1_init, η, maxiter)
ulp(w1) = log_joint(w1, params...)
w1_seq = gradient_method_1dim(ulp, w1_init, η, maxiter)
w1_seq
end
# 対数同時分布
log_joint(w1, X, Y, w2, σ, μ1, σ1) = sum(logpdf.(Normal.(w1*X .+ w2, σ), Y)) + logpdf(Normal(μ1, σ1), w1)
params = (X_obs, Y_obs, w2, σ, μ1, σ1)
# 最適化の実施
w1_seq = inference_wrapper_gd_1dim(log_joint, params, w1_init, η, maxiter)
# 勾配法の過程を可視化
fig, ax = subplots(figsize=(8,4))
ax.plot(w1_seq)
set_options(ax, "iteration", "w1", "w1 sequence")
false
# 最適化パラメータ
w1_init = 0.0
maxiter = 100
η = 0.01
# グローバル変数を参照する場合
@time gradient_method_1dim(ulp, w1_init, η, maxiter)
@time inference_wrapper_gd_1dim(log_joint, params, w1_init, η, maxiter)
0.000263 seconds (1.98 k allocations: 76.672 KiB) 0.000042 seconds (796 allocations: 35.094 KiB)
100-element Vector{Float64}:
0.0
0.109
0.18528909999999998
0.23868384109
0.276054820378891
0.3022107687831858
0.3205173170713517
0.33333007021823907
0.3422977161457455
0.34857417153040726
0.352967062654132
0.356041647151627
0.35819354884142374
⋮
0.3632122625791319
0.3632122625791344
0.3632122625791362
0.3632122625791374
0.3632122625791383
0.3632122625791389
0.3632122625791393
0.36321226257913963
0.36321226257913986
0.36321226257913997
0.3632122625791401
0.36321226257914013
# 近似分布用の平均を求める
μ_approx = w1_seq[end]
# 分布を表示する範囲
w1s = range(-5, 5, length=100)
fig, axes = subplots(1, 2, figsize=(8,4))
# 非正規化対数事後分布と最適値
axes[1].plot(w1s, ulp.(w1s))
axes[1].plot(μ_approx, ulp(μ_approx), "rx", label="optimal")
set_options(axes[1], "w1", "log density (unnormalized)", "unnormalized log posterior distribution")
# 非正規化事後分布と最適値
axes[2].plot(w1s, exp.(ulp.(w1s)))
axes[2].plot(μ_approx, exp.(ulp(μ_approx)), "rx", label="optimal")
set_options(axes[2], "w1", "density (unnormalized)", "unnormalized posterior distribution")
tight_layout()
# 2階微分から分散を求める
grad(x) = ForwardDiff.derivative(ulp, x)
hessian(x) = ForwardDiff.derivative(grad, x)
σ_approx = sqrt(inv(-hessian(μ_approx)))
0.1825437644092281
# 可視化する範囲
w1s = range(-0.5, 1.0, length=100)
# 非正規化事後分布の可視化
fig, axes = subplots(1, 2, figsize=(8,4))
axes[1].plot(w1s, exp.(ulp.(w1s)))
set_options(axes[1], "w1", "density (unnormalized)", "unnormalized posterior distributions")
# 得られた近似分布の可視化
axes[2].plot(w1s, pdf.(Normal.(μ_approx, σ_approx), w1s))
set_options(axes[2], "w1", "density", "approximate distribution")
tight_layout()
# 事前分布の設定
σ = 1.0
μ1 = 0.0
μ2 = 0.0
σ1 = 10.0
σ2 = 10.0
# 対数同時分布
log_joint(w, X, Y, σ, μ1, σ1, μ2, σ2) =
sum(logpdf.(Normal.(w[1]*X .+ w[2], σ), Y)) +
logpdf(Normal(μ1, σ1), w[1]) +
logpdf(Normal(μ2, σ2), w[2])
# 非正規化対数事後分布
params = (X_obs, Y_obs, σ, μ1, σ1, μ2, σ2)
ulp(w) = log_joint(w, params...)
# 分布を可視化する範囲
w1s = range(-5, 5, length=100)
w2s = range(-5, 5, length=100)
fig, axes = subplots(1, 2, figsize=(12,4))
# 非正規化対数事後分布の可視化
cs = axes[1].contour(w1s, w2s, [ulp([w1, w2]) for w1 in w1s, w2 in w2s]', cmap="jet")
axes[1].clabel(cs, inline=true)
set_options(axes[1], "w1", "w2", "unnormalized log posterior")
# 非正規化事後分布の可視化
cs = axes[2].contour(w1s, w2s, [exp(ulp([w1, w2])) for w1 in w1s, w2 in w2s]', cmap="jet")
axes[2].clabel(cs, inline=true)
set_options(axes[2], "w1", "w2", "unnormalized posterior")
tight_layout()
# 多次元の勾配法
function gradient_method(funct, x_init, η, maxiter)
x_seq = Array{typeof(x_init[1]), 2}(undef, length(x_init), maxiter)
g(x) = ForwardDiff.gradient(funct, x)
x_seq[:, 1] = x_init
for i in 2:maxiter
x_seq[:,i] = x_seq[:, i-1] + η*g(x_seq[:, i-1])
end
x_seq
end
function inference_wrapper_gd(log_joint, params, w_init, η, maxiter)
ulp(w) = log_joint(w, params...)
w_seq = gradient_method(ulp, w_init, η, maxiter)
w_seq
end
# 最適化パラメータ
w_init = [0.0, 0.0]
maxiter = 1000
η = 0.01
# 最適化の実行
w_seq = inference_wrapper_gd(log_joint, params, w_init, η, maxiter)
# 勾配法の過程を可視化
fig, axes =subplots(2, 1, figsize=(8,8))
axes[1].plot(w_seq[1, :])
set_options(axes[1], "iteration", "w1", "w1 sequence")
axes[2].plot(w_seq[2,:])
set_options(axes[2], "iteration", "w2", "w2, sequence")
tight_layout()
# 平均
μ_approx = w_seq[:, end]
# 共分散
hessian(w) = ForwardDiff.hessian(ulp, w)
Σ_approx = inv(-hessian(μ_approx))
fog, axes =subplots(1, 2, figsize=(8,4))
# 非正規化事後分布の可視化
cs = axes[1].contour(w1s, w2s, [exp(ulp([w1, w2])) for w1 in w1s, w2 in w2s]', cmap="jet")
axes[1].clabel(cs, inline=true)
set_options(axes[1], "w1", "w2", "unnormalized posterior")
# 近似正規分布の可視化
cx = axes[2].contour(w1s, w2s, [pdf(MvNormal(μ_approx, Σ_approx), [w1, w2]) for w1 in w1s, w2 in w2s]', cmap="jet")
axes[2].clabel(cs, inline=true)
set_options(axes[2], "w1", "w2", "approximate posterior")
tight_layout()
# 長方形の幅
Δ1 = w1s[2] - w1s[1]
Δ2 = w2s[2] - w2s[1]
# 積分近似による予測分布
p_predictive2(x, y) = sum([pdf(Normal(w1*x + w2, σ), y) *
pdf(MvNormal(μ_approx, Σ_approx), [w1, w2]) *
Δ1 * Δ2 for w1 in w1s, w2 in w2s])
# 描画範囲
xs = range(-10, 10, length=100)
ys = range(-5, 5, length=100)
# 密度の計算
density_y = p_predictive2.(xs, ys')
fig, axes = subplots(1, 2, sharey=true, figsize=(8,4))
# 等高線図
cs = axes[1].contour(xs, ys, density_y', cmap="jet")
axes[1].clabel(cs, inline=true)
axes[1].scatter(X_obs, Y_obs)
set_options(axes[1], "x", "y", "predictive distribution (contour)")
# カラーメッシュ
xgrid = repeat(xs', length(ys), 1)
ygrid = repeat(ys, 1, length(xs))
axes[2].pcolormesh(xgrid, ygrid, density_y', cmap="jet", shading="auto")
axes[2].plot(X_obs, Y_obs, "ko", label="data")
set_options(axes[2], "x", "y", "predictive distribution (colored)")
tight_layout()
# 入力データ集合
X_obs = [-2, 1, 2]
# 出力データ
Y_obs = Bool.([0, 1, 1])
# 散布図で可視化
fig, ax = subplots()
ax.scatter(X_obs, Y_obs)
ax.set_xlim([-3, 3])
set_options(ax, "x", "y", "data (N = $(length(X_obs)))")
false
# シグモイド関数
sig(x) = 1/(1 + exp(-x))
# 事前分布の設定
μ1 = 0
μ2 = 0
σ1 = 10.0
σ2 = 10.0
# 対数同時分布
log_joint(w, X, Y, μ1, σ1, μ2, σ2) =
logpdf(Normal(μ1, σ1), w[1]) +
logpdf(Normal(μ2, σ2), w[2]) +
sum(logpdf.(Bernoulli.(sig.(w[1]*X .+ w[2])), Y))
params = (X_obs, Y_obs, μ1, σ1, μ2, σ2)
# 非正規化対数事後分布
ulp(w) = log_joint(w, params...)
ulp (generic function with 1 method)
# 最適化パラメータ
w_init = [0.0, 0.0]
maxiter = 2000
η = 0.1
# 最適化の実行
w_seq = inference_wrapper_gd(log_joint, params, w_init, η, maxiter)
# 勾配法の過程を可視化
fig, axes = subplots(2, 1, figsize=(8,8))
axes[1].plot(w_seq[1,:])
set_options(axes[1], "iteration", "w1", "w1 sequence")
axes[2].plot(w_seq[2,:])
set_options(axes[2], "iteration", "w2", "w2 sequence")
tight_layout()
# 平均
μ_approx = w_seq[:, end]
# 共分散
hessian(x) = ForwardDiff.hessian(ulp, x)
Σ_approx = inv(-hessian(μ_approx))
2×2 Matrix{Float64}:
18.6137 -2.69261
-2.69261 30.6517
w1s = range(-10, 30, length=100)
w2s = range(-20, 20, length=100)
fig, axes = subplots(1, 2, figsize=(8,4))
cs = axes[1].contour(w1s, w2s, [exp(ulp([w1, w2])) + eps() for w1 in w1s, w2 in w2s]', cmap="jet")
axes[1].clabel(cs, inline=true)
set_options(axes[1], "w1", "w2", "unnormalized posterior")
cs = axes[2].contour(w1s, w2s, [pdf(MvNormal(μ_approx, Σ_approx), [w1, w2]) for w1 in w1s, w2 in w2s]', cmap="jet")
axes[2].clabel(cs, inline=true)
set_options(axes[2], "w1", "w2", "approximate posterior")
tight_layout()
# 近似分布から候補のサンプルを１００個抽出
W = rand(MvNormal(μ_approx, Σ_approx), 100)
fig, ax = subplots()
for i in 1:size(W, 2)
w1, w2 = W[:, i]
funct(x) = sig(w1*x + w2)
ax.plot(xs, funct.(xs), "g", alpha=0.1)
end
ax.scatter(X_obs, Y_obs)
ax.set_xlim(extrema(xs))
set_options(ax, "x", "y", "prediction sample from approximate posterior")
false
# 長方形の幅
Δ1 = w1s[2] - w1s[1]
Δ2 = w2s[2] - w2s[1]
# 積分近似による予測分布
p_predictive2(x, y) = sum([pdf(Bernoulli(sig(w1*x + w2)), y) *
pdf(MvNormal(μ_approx, Σ_approx), [w1, w2]) *
Δ1 * Δ2 for w1 in w1s, w2 in w2s])
# 関数を可視化する範囲
xs = range(-3, 3, length=50)
fig, ax = subplots()
ax.scatter(X_obs, Y_obs, label="data")
ax.plot(xs, p_predictive2.(xs, 1), label="probability")
ax.set_xlim([-3, 3])
set_options(ax, "x", "y", "predection", legend=true)
PyObject <matplotlib.legend.Legend object at 0x7fe03ef0ae20>
# ガウス提案分布によるメトロポリスヘイスティング
function GaussianMH(log_p_tilde, μ0; maxiter::Int=100_000, σ::Float64=1.0)
# サンプルを格納する配列
D = length(μ0)
μ_samples = Array{typeof(μ0[1]), 2}(undef, D, maxiter)
# 初期サンプル
μ_samples[:, 1] = μ0
# 受容されたサンプルの数
num_accepted = 1
for i in 2:maxiter
# 提案分布q（多変量正規分布）に従い、次のサンプルの候補を抽出
μ_tmp = rand(MvNormal(μ_samples[:, i-1], σ*eye(D)))
# 比率r（の対数）を計算
log_r = (log_p_tilde(μ_tmp) + logpdf(MvNormal(μ_tmp, σ*eye(D)), μ_samples[:, i-1])) -
(log_p_tilde(μ_samples[:, i-1]) + logpdf(MvNormal(μ_samples[:, i-1], σ*eye(D)), μ_tmp))
# 確率rでサンプルを受容する
is_accepted = min(1, exp(log_r)) > rand()
new_sample = is_accepted ? μ_tmp : μ_samples[:, i-1]
# 新しいサンプルを格納
μ_samples[:, i] = new_sample
# 受容された場合、合計をプラスする
num_accepted += is_accepted
end
μ_samples, num_accepted
end
GaussianMH (generic function with 1 method)
# ハミルトニアンモンテカルロ法
function HMC(log_p_tilde, μ0; maxiter::Int=100_000, L::Int=100, ϵ::Float64=1e-1)
# leapflogによる値の更新
function leapflog(grad, p_in, μ_in, L, ϵ)
μ = μ_in
p = p_in + 0.5*ϵ*grad(μ)
for l in 1:L-1
μ += ϵ*p
p += ϵ*grad(μ)
end
μ += ϵ*p
p += 0.5*ϵ*grad(μ)
p, μ
end
# 非正規化対数事後分布の勾配関数の計算
grad(μ) = ForwardDiff.gradient(log_p_tilde, μ)
# サンプルを格納する配列
D = length(μ0)
μ_samples = Array{typeof(μ0[1]), 2}(undef, D, maxiter)
# 初期サンプル
μ_samples[:, 1] = μ0
# 受容されたサンプルの数
num_accepted = 1
for i in 2:maxiter
# 運動量pの生成
p_in = randn(size(μ0))
# リープフロッグ法
p_out, μ_out = leapflog(grad, p_in, μ_samples[:, i-1], L, ϵ)
# 比率r（の対数）を計算
μ_in = μ_samples[:, i-1]
log_r = (log_p_tilde(μ_out) + logpdf(MvNormal(zeros(D), eye(D)), vec(p_out))) -
(log_p_tilde(μ_in) + logpdf(MvNormal(zeros(D), eye(D)), vec(p_in)))
# 確率rでサンプルを受容する
is_accepted = min(1, exp(log_r)) > rand()
new_sample = is_accepted ? μ_out : μ_in
# 新しいサンプルを受容する
μ_samples[:, i] = new_sample
# 受容された場合、合計をプラスする
num_accepted += is_accepted
end
μ_samples, num_accepted
end
HMC (generic function with 1 method)
function inference_wrapper_GMH(log_joint, params, w_init; maxiter::Int=100_000, σ::Float64=1.0)
ulp(w) = log_joint(w, params...)
GaussianMH(ulp, w_init; maxiter=maxiter, σ=σ)
end
function inference_wrapper_HMC(log_joint, params, w_init; maxiter::Int=100_000, L::Int=100, ϵ::Float64=1e-1)
ulp(w) = log_joint(w, params...)
HMC(ulp, w_init; maxiter=maxiter, L=L, ϵ=ϵ)
end
inference_wrapper_HMC (generic function with 1 method)
# 入力データセット
X_obs = [-2, 1, 5]
# 出力データセット
Y_obs = [-2.2, -1.0, 1.5]
# 散布図で可視化
fig, ax = subplots()
ax.scatter(X_obs, Y_obs)
set_options(ax, "x", "y", "data (N = $(length(X_obs)))")
false
# 対数同時分布
log_joint(w, X, Y, σ, μ1, σ1, μ2, σ2) =
sum(logpdf.(Normal.(w[1]*X .+ w[2], σ), Y)) +
logpdf(Normal(μ1, σ1), w[1]) +
logpdf(Normal(μ2, σ2), w[2])
params = (X_obs, Y_obs, σ, μ1, σ1, μ2, σ2)
# 非正規化対数事後分布
ulp(w) = log_joint(w, params...)
ulp (generic function with 1 method)
# 初期値
w_init = randn(2)
# サンプリング
maxiter = 300
param_posterior_GMH, num_accepted_GMH =
inference_wrapper_GMH(log_joint, params, w_init, maxiter=maxiter, σ=1.0)
param_posterior_HMC, num_accepted_HMC =
inference_wrapper_HMC(log_joint, params, w_init, maxiter=maxiter, L=10, ϵ=1e-1)
# サンプリングの過程を可視化（Gaussian MH）
fig, axes = subplots(2, 1, figsize=(8,4))
axes[1].plot(param_posterior_GMH[1, :])
set_options(axes[1], "iteration", "w1", "w1 sequence (GMH)")
axes[2].plot(param_posterior_GMH[2,:])
set_options(axes[2], "iteration", "w2", "w2 sequence (GMH)")
tight_layout()
println("acceptance rate (GMH) = $(num_accepted_GMH/maxiter)")
# サンプリングの過程を可視化（HMC）
fig, axes = subplots(2, 1, figsize=(8,4))
axes[1].plot(param_posterior_HMC[1, :])
set_options(axes[1], "iteration", "w1", "w1 sequence (HMC)")
axes[2].plot(param_posterior_HMC[2,:])
set_options(axes[2], "iteration", "w2", "w2 sequence (HMC)")
tight_layout()
println("acceptance rate (HMC) = $(num_accepted_HMC/maxiter)")
acceptance rate (GMH) = 0.16 acceptance rate (HMC) = 0.9933333333333333
# 事後分布を可視化する範囲
w1s = range(-5, 5, length=100)
w2s = range(-5, 5, length=100)
fig, axes =subplots(1, 3, sharex=true, sharey=true, figsize=(12,4))
cs = axes[1].contour(w1s, w2s, [exp(ulp([w1, w2])) + eps() for w1 in w1s, w2 in w2s]', cmap="jet")
# 非正規化事後分布
axes[1].clabel(cs, inline=true)
set_options(axes[1], "w1", "w2", "unnormalized posterior")
# Gaussian MHで得られた事後分布からのサンプル
axes[2].scatter(param_posterior_GMH[1, :], param_posterior_GMH[2, :], alpha=10/maxiter)
set_options(axes[2], "w1", "w2", "samples from posterior (GMH)")
# HMCで得られた事後分布からのサンプル
axes[3].scatter(param_posterior_HMC[1, :], param_posterior_HMC[2, :], alpha=10/maxiter)
set_options(axes[3], "w1", "w2", "samples from posterior (HMC)")
false
xs = range(-10, 10, length=100)
fig, axes = subplots(1, 2, figsize=(8,4))
# Gaussian MH
for i in 1:size(param_posterior_GMH, 2)
w1, w2 = param_posterior_GMH[:, i]
funct(x) = w1*x + w2
axes[1].plot(xs, funct.(xs), "g", alpha=10/maxiter)
end
axes[1].scatter(X_obs, Y_obs)
axes[1].set_xlim(extrema(xs))
set_options(axes[1], "x", "y", "predictive distribution (GMH)")
# HMC
for i in 1:size(param_posterior_HMC, 2)
w1, w2 = param_posterior_HMC[:, i]
funct(x) = w1*x + w2
axes[2].plot(xs, funct.(xs), "g", alpha=10/maxiter)
end
axes[2].scatter(X_obs, Y_obs)
axes[2].set_xlim(extrema(xs))
set_options(axes[2], "x", "y", "predictive distribution (HMC)")
false
# 入力データセット
X_obs = [-2, 1, 2]
# 出力データセット
Y_obs = Bool.([0, 1, 1])
# 散布図で可視化
fig, ax = subplots()
ax.scatter(X_obs, Y_obs)
set_options(ax, "x", "y", "data = (N = $(length(X_obs)))")
false
# シグモイド関数
sig(x) = 1/(1 + exp(-x))
# 事前分布の設定
μ1 = 0
μ2 = 0
σ1 = 10.0
σ2 = 10.0
# 対数同時分布
log_joint(w, X, Y, μ1, σ1, μ2, σ2) =
logpdf(Normal(μ1, σ1), w[1]) +
logpdf(Normal(μ2, σ2), w[2]) +
sum(logpdf.(Bernoulli.(sig.(w[1]*X .+ w[2])), Y))
params = (X_obs, Y_obs, μ1, σ1, μ2, σ2)
# 非正規化対数事後分布
ulp(w) = log_joint(w, params...)
ulp (generic function with 1 method)
# 初期値
w_init = randn(2)
# サンプリング
maxiter = 300
param_posterior_GMH, num_accepted_GMH =
inference_wrapper_GMH(log_joint, params, w_init, maxiter=maxiter, σ=1.0)
param_posterior_HMC, num_accepted_HMC =
inference_wrapper_HMC(log_joint, params, w_init, maxiter=maxiter, L=10, ϵ=1e-1)
# サンプリングの過程を可視化（Gaussian MH）
fig, axes = subplots(2, 1, figsize=(8,4))
axes[1].plot(param_posterior_GMH[1, :])
set_options(axes[1], "iteration", "w1", "w1 sequence (GMH)")
axes[2].plot(param_posterior_GMH[2,:])
set_options(axes[2], "iteration", "w2", "w2 sequence (GMH)")
tight_layout()
println("acceptance rate (GMH) = $(num_accepted_GMH/maxiter)")
# サンプリングの過程を可視化（HMC）
fig, axes = subplots(2, 1, figsize=(8,4))
axes[1].plot(param_posterior_HMC[1, :])
set_options(axes[1], "iteration", "w1", "w1 sequence (HMC)")
axes[2].plot(param_posterior_HMC[2,:])
set_options(axes[2], "iteration", "w2", "w2 sequence (HMC)")
tight_layout()
println("acceptance rate (HMC) = $(num_accepted_HMC/maxiter)")
acceptance rate (GMH) = 0.9033333333333333 acceptance rate (HMC) = 1.0
# 事後分布を可視化する範囲
w1s = range(-10, 30, length=100)
w2s = range(-20, 20, length=100)
fig, axes =subplots(1, 3, sharex=true, sharey=true, figsize=(12,4))
cs = axes[1].contour(w1s, w2s, [exp(ulp([w1, w2])) + eps() for w1 in w1s, w2 in w2s]', cmap="jet")
# 非正規化事後分布
axes[1].clabel(cs, inline=true)
set_options(axes[1], "w1", "w2", "unnormalized posterior")
# Gaussian MHで得られた事後分布からのサンプル
axes[2].scatter(param_posterior_GMH[1, :], param_posterior_GMH[2, :], alpha=10/maxiter)
set_options(axes[2], "w1", "w2", "samples from posterior (GMH)")
# HMCで得られた事後分布からのサンプル
axes[3].scatter(param_posterior_HMC[1, :], param_posterior_HMC[2, :], alpha=10/maxiter)
set_options(axes[3], "w1", "w2", "samples from posterior (HMC)")
false
param_posterior_GMH, num_accepted_GMH =
GaussianMH(ulp, w_init, maxiter=maxiter, σ=10.0)
param_posterior_HMC, num_accepted_HMC =
HMC(ulp, w_init, maxiter=maxiter, L=10, ϵ=1e-0)
# サンプリングの過程を可視化（Gaussian MH）
fig, axes = subplots(2, 1, figsize=(8,4))
axes[1].plot(param_posterior_GMH[1, :])
set_options(axes[1], "iteration", "w1", "w1 sequence (GMH)")
axes[2].plot(param_posterior_GMH[2,:])
set_options(axes[2], "iteration", "w2", "w2 sequence (GMH)")
tight_layout()
println("acceptance rate (GMH) = $(num_accepted_GMH/maxiter)")
# サンプリングの過程を可視化（HMC）
fig, axes = subplots(2, 1, figsize=(8,4))
axes[1].plot(param_posterior_HMC[1, :])
set_options(axes[1], "iteration", "w1", "w1 sequence (HMC)")
axes[2].plot(param_posterior_HMC[2,:])
set_options(axes[2], "iteration", "w2", "w2 sequence (HMC)")
tight_layout()
println("acceptance rate (HMC) = $(num_accepted_HMC/maxiter)")
acceptance rate (GMH) = 0.6466666666666666 acceptance rate (HMC) = 0.9766666666666667
# 事後分布を可視化する範囲
w1s = range(-10, 30, length=100)
w2s = range(-20, 20, length=100)
fig, axes =subplots(1, 3, sharex=true, sharey=true, figsize=(12,4))
cs = axes[1].contour(w1s, w2s, [exp(ulp([w1, w2])) + eps() for w1 in w1s, w2 in w2s]', cmap="jet")
# 非正規化事後分布
axes[1].clabel(cs, inline=true)
set_options(axes[1], "w1", "w2", "unnormalized posterior")
# Gaussian MHで得られた事後分布からのサンプル
axes[2].scatter(param_posterior_GMH[1, :], param_posterior_GMH[2, :], alpha=10/maxiter)
set_options(axes[2], "w1", "w2", "samples from posterior (GMH)")
# HMCで得られた事後分布からのサンプル
axes[3].scatter(param_posterior_HMC[1, :], param_posterior_HMC[2, :], alpha=10/maxiter)
set_options(axes[3], "w1", "w2", "samples from posterior (HMC)")
false
# 関数を可視化する範囲
xs = range(-3, 3, length=100)
fig, axes = subplots(2, 2, figsize=(12,8))
# Gaussian MHによるサンプル
fs =[]
for i in 1:size(param_posterior_GMH, 2)
w1, w2 = param_posterior_GMH[:, i]
funct(x) = sig(w1*x + w2)
push!(fs, funct.(xs))
axes[1,1].plot(xs, funct.(xs), "g", alpha=10/maxiter)
end
axes[1,1].scatter(X_obs, Y_obs)
axes[1,1].set_xlim(extrema(xs))
set_options(axes[1,1], "x", "y", "predictive distribution (GMH)")
# GMHによる予測の平均
axes[1,2].plot(xs, mean(fs), label="prediction")
axes[1,2].scatter(X_obs, Y_obs, label="data")
set_options(axes[1,2], "x", "y (prob.)", "prediction (GMH)", legend=true)
# HMCによるサンプル
fs =[]
for i in 1:size(param_posterior_HMC, 2)
w1, w2 = param_posterior_HMC[:, i]
funct(x) = sig(w1*x + w2)
push!(fs, funct.(xs))
axes[2,1].plot(xs, funct.(xs), "g", alpha=10/maxiter)
end
axes[2,1].scatter(X_obs, Y_obs)
axes[2,1].set_xlim(extrema(xs))
set_options(axes[2,1], "x", "y", "predictive distribution (HMC)")
# HMCによる予測の平均
axes[2,2].plot(xs, mean(fs), label="prediction")
axes[2,2].scatter(X_obs, Y_obs, label="data")
set_options(axes[2,2], "x", "y (prob.)", "prediction (HMC)", legend=true)
PyObject <matplotlib.legend.Legend object at 0x7fe055642130>
# 入力データセット
X_obs = [-2, -1.5, 0.5, 0.7, 1.2]
# 出力データセット
Y_obs = [0, 0, 2, 1, 8]
# 散布図で可視化
fig, ax = subplots()
ax.scatter(X_obs, Y_obs)
ax.set_xlim([-2, 2])
set_options(ax, "x", "y", "data = (N = $(length(X_obs)))")
false
# 対数同時分布
log_joint(w, X, Y) =
sum(logpdf.(Poisson.(exp.(w[1]*X .+ w[2])), Y)) +
logpdf(Normal(μ1, σ1), w1) +
logpdf(Normal(μ2, σ2), w2)
params = (X_obs, Y_obs)
# 非正規化対数事後分布
ulp(w) = log_joint(w, params...)
# 初期値
w_init = randn(2)
# サンプルサイズ
maxiter = 300
# ハミルトニアンモンテカルロ法によるサンプリング
param_posterior, num_accepted = inference_wrapper_HMC(log_joint, params, w_init, maxiter=maxiter)
# トレースプロット
fig, axes = subplots(2, 1, figsize=(8,4))
axes[1].plot(param_posterior[1, :])
set_options(axes[1], "iteration", "w1", "w1 sequence")
axes[2].plot(param_posterior[2, :])
set_options(axes[2], "iteration", "w2", "w2 sequence")
tight_layout()
println("acceptance rate = $(num_accepted/maxiter)")
# HMCで得られた事後分布からのサンプル
fig, ax = subplots()
ax.scatter(param_posterior[1, :], param_posterior[2, :], alpha=0.1)
set_options(ax, "w1", "w2", "sample from posterior")
tight_layout()
acceptance rate = 0.97
# 関数を可視化する範囲
xs = range(-2, 2, length=100)
fig, ax = subplots()
for i in 1:size(param_posterior, 2)
w1, w2 = param_posterior[:, i]
# 指数関数を可視化
funct(x) = exp(w1*x + w2)
ax.plot(xs, funct.(xs), "g", alpha=0.1)
end
ax.plot(X_obs, Y_obs, "ko")
ax.set_xlim(extrema(xs))
ax.set_ylim([-1, 15])
set_options(ax, "x", "y", "predictive distribution")
tight_layout()
# 学習データ
X_obs = []
Y_obs = []
push!(X_obs, [0.3, 0.4])
push!(Y_obs, [4.0, 3.7])
push!(X_obs, [0.2, 0.4, 0.9])
push!(Y_obs, [6.0, 7.2, 9.4])
push!(X_obs, [0.6, 0.8, 0.9])
push!(Y_obs, [6.0, 6.9, 7.8])
# 散布図の可視化
fig, ax = subplots()
ax.plot(X_obs[1], Y_obs[1], "or", label="class1")
ax.plot(X_obs[2], Y_obs[2], "xg", label="class2")
ax.plot(X_obs[3], Y_obs[3], "^b", label="class3")
set_options(ax, "x", "y", "data", legend=true)
PyObject <matplotlib.legend.Legend object at 0x7fe03ea39e80>
function linear_fit(Y, X)
N = length(Y)
w1 = sum((Y .- mean(Y)) .* X) / sum((X .- mean(X)) .* X)
w2 = mean(Y) - w1*mean(X)
w1, w2
end
# まとめて回帰
w1, w2 = linear_fit(vcat(Y_obs...), vcat(X_obs...))
# 個別に回帰
w1s = []
w2s = []
for i in 1:3
w1_tmp, w2_tmp = linear_fit(Y_obs[i], X_obs[i])
push!(w1s, w1_tmp)
push!(w2s, w2_tmp)
end
# 関数を可視化する範囲
xs = range(0, 1, length=100)
fig, axes = subplots(1, 2, figsize=(8,4))
# 単一の回帰
axes[1].plot(xs, w1.*xs .+ w2, "-k")
# 個別の回帰
axes[2].plot(xs, w1s[1].*xs .+ w2s[1], "-r")
axes[2].plot(xs, w1s[2].*xs .+ w2s[2], "-g")
axes[2].plot(xs, w1s[3].*xs .+ w2s[3], "-b")
# データの可視化
for ax in axes
ax.plot(X_obs[1], Y_obs[1], "or", label="class1")
ax.plot(X_obs[2], Y_obs[2], "xg", label="class2")
ax.plot(X_obs[3], Y_obs[3], "^b", label="class3")
end
set_options(axes[1], "x", "y", "(a) single regression", legend=true)
set_options(axes[2], "x", "y", "(b) multiple regressions", legend=true)
tight_layout()
# 対数同時分布の設計
@views hyper_prior(w) = logpdf(Normal(0, 10.0), w[1]) + logpdf(Normal(0, 10.0), w[2])
@views prior(w) = sum(logpdf.(Normal.(w[1], 1.0), w[3:5])) + sum(logpdf.(Normal.(w[2], 1.0), w[6:8]))
@views log_likelihood(Y, X, w) = sum([sum(logpdf.(Normal.(Y[i], 1.0), w[2+i].*X[i] .+ w[2+i+3])) for i in 1:3])
log_joint(w, X, Y) = hyper_prior(w) + prior(w) + log_likelihood(Y, X, w)
params = (X_obs, Y_obs)
ulp(w) = log_joint(w, params...)
ulp (generic function with 1 method)
# HMCによるサンプリング
maxiter = 1000
w_init = randn(8)
param_posterior, num_accepted = inference_wrapper_HMC(log_joint, params, w_init, maxiter=maxiter, L=100, ϵ=1e-2)
# 予測分布の可視化
fig, axes = subplots(1, 3, sharey=true, figsize=(12, 4))
for i in 1:3
for j in 1:size(param_posterior, 2)
w1, w2 = param_posterior[[2+i, 2+i+3], j]
axes[i].plot(xs, w1.*xs .+ w2, "c-", alpha=0.01)
end
set_options(axes[i], "x", "y", "class $(i)")
end
axes[1].plot(X_obs[1], Y_obs[1], "or")
axes[2].plot(X_obs[2], Y_obs[2], "xg")
axes[3].plot(X_obs[3], Y_obs[3], "^b")
tight_layout()
X_obs = []
Y_obs = []
push!(X_obs, [0.1, 0.3, 0.4, 0.5, 0.6, 0.9])
push!(Y_obs, [4.0, 4.0, 3.7, 3.8, 3.9, 3.7])
push!(X_obs, [0.2, 0.4, 0.9])
push!(Y_obs, [6.0, 7.2, 9.4])
push!(X_obs, [0.6, 0.8, 0.9])
push!(Y_obs, [6.0, 6.9, 7.8])
# 対数同時分布の設計
@views hyper_prior(w) = logpdf(Normal(0, 10.0), w[1]) + logpdf(Normal(0, 10.0), w[2])
@views prior(w) = sum(logpdf.(Normal.(w[1], 1.0), w[3:5])) + sum(logpdf.(Normal.(w[2], 1.0), w[6:8]))
@views log_likelihood(Y, X, w) = sum([sum(logpdf.(Normal.(Y[i], 1.0), w[2+i].*X[i] .+ w[2+i+3])) for i in 1:3])
log_joint(w, X, Y) = hyper_prior(w) + prior(w) + log_likelihood(Y, X, w)
params = (X_obs, Y_obs)
ulp(w) = log_joint(w, params...)
# HMCによるサンプリング
maxiter = 1000
w_init = randn(8)
param_posterior, num_accepted = inference_wrapper_HMC(log_joint, params, w_init, maxiter=maxiter, L=100, ϵ=1e-2)
# 予測分布の可視化
fig, axes = subplots(1, 3, sharey=true, figsize=(12, 4))
for i in 1:3
for j in 1:size(param_posterior, 2)
w1, w2 = param_posterior[[2+i, 2+i+3], j]
axes[i].plot(xs, w1.*xs .+ w2, "c-", alpha=0.01)
end
set_options(axes[i], "x", "y", "class $(i)")
end
axes[1].plot(X_obs[1], Y_obs[1], "or")
axes[2].plot(X_obs[2], Y_obs[2], "xg")
axes[3].plot(X_obs[3], Y_obs[3], "^b")
tight_layout()
# 系列の長さ
N = 20
# 観測データの次元
D = 2
# 系列データ（#=, =#は複数業をつなげるために挿入）
Y_obs =
[1.9 0.2 0.1 1.4 0.3 1.3 1.6 1.5 1.6 2.4 #=
=# 1.7 3.6 2.8 1.6 3.0 2.8 5.1 5.2 6.0 6.4;
0.1 0.2 0.9 1.5 4.0 5.0 6.3 5.8 6.4 7.5 #=
=# 6.7 7.6 8.7 8.2 8.5 9.6 8.4 8.4 8.4 9.0]
# ２次元に系列データを可視化
fig, ax = subplots(figsize=(6,6))
ax.plot(Y_obs[1, :], Y_obs[2, :], "kx--")
ax.text(Y_obs[1, 1], Y_obs[2, 1], "start", color="r")
ax.text(Y_obs[1, end], Y_obs[2, end], "goal", color="r")
set_options(ax, "y1", "y2", "2dim data")
false
# 初期状態に与えるノイズ量
σ1 = 100.0
# 状態の遷移に仮定するノイズ量
σ_x = 1.0
# 観測に仮定するノイズ量
σ_y = 1.0
# 状態の遷移系列に関する対数密度
@views transition(X, σ1, σ_x, D, N) =
logpdf(MvNormal(zeros(D), σ1 * eye(D)), X[:, 1]) +
sum([logpdf(MvNormal(X[:, n-1], σ_x * eye(D)), X[:, n]) for n in 2:N])
# 観測データに関する対数密度
@views observation(X, Y, σ_y, D, N) =
sum([logpdf(MvNormal(X[:, n], σ_y * eye(D)), Y[:, n]) for n in 1:N])
# 対数同時分布
log_joint_tmp(X, Y, σ1, σ_x, σ_y, D, N) =
transition(X, σ1, σ_x, D, N) +
observation(X, Y, σ_y, D, N)
# DN次元ベクトルを入力とする関数にする
log_joint(X_vec, Y, σ1, σ_x, σ_y, D, N) =
log_joint_tmp(reshape(X_vec, D, N), Y, σ1, σ_x, σ_y, D, N)
params = (Y_obs, σ1, σ_x, σ_y, D, N)
# 非正規化対数事後分布
ulp(X_vec) = log_joint(X_vec, params...)
ulp (generic function with 1 method)
# 初期値
X_init = randn(D * N)
# サンプルサイズ
maxiter = 1000
# HMCの実行
samples, num_accepted = inference_wrapper_HMC(log_joint, params, X_init, maxiter=maxiter)
println("acceptance rate = $(num_accepted/maxiter)")
acceptance rate = 0.995
fig, ax = subplots(figsize=(6,6))
for i in 1:maxiter
# i番目のサンプルの状態変数
X = reshape(samples[:, i], D, N)
ax.plot(X[1, :], X[2, :], "go--", alpha=10.0/maxiter)
end
# 観測データの可視化
ax.plot(Y_obs[1, :], Y_obs[2, :], "kx--", label="observation (Y)")
# 状態変数の平均
mean_trace = reshape(mean(samples, dims=2), D, N)
ax.plot(mean_trace[1, :], mean_trace[2, :], "ro--", label="estimated position (X)")
set_options(ax, "y1", "y2", "2dim data", legend=true)
PyObject <matplotlib.legend.Legend object at 0x7fe041bbb8b0>
# 観測データ数
N = 20
# 入力値
Z_obs = [10, 10, 10, 10, 10, 10, 10, 10, 10, 15, 15, 15, 15, 15, 15, 15, 8, 8, 8, 8]
# 出力値
Y_obs = [67, 64, 60, 60, 57, 54, 51, 51, 49, 63, 62, 62, 58, 57, 53, 51, 24, 22, 23, 19]
# データの可視化
fig, ax = subplots()
ax.scatter(Z_obs, Y_obs)
set_options(ax, "z", "y", "data (scatter)")
false
fig, axes = subplots(2, 1, figsize=(8,6))
# 出力値
axes[1].plot(Y_obs)
set_options(axes[1], "time", "y", "time series (Y_obs)")
# 入力値
axes[2].plot(Z_obs)
set_options(axes[2], "time", "z", "time series (Z_obs)")
false
# 初期状態に与えるノイズ量
σ1 = 10.0
# 状態の遷移に仮定するノイズ量
σ_x = 1.0
# 観測に仮定するノイズ量
σ_y = 0.5
# パラメータに仮定するノイズ量
σ_w = 100.0
# パラメータの事前分布
prior(w, σ_w) = logpdf(MvNormal(zeros(2), σ_w*eye(2)), w)
# 状態の遷移系列に関する対数密度
@views transition(X, σ1, σ_x, N) =
logpdf(Normal(0, σ1), X[1]) +
sum(logpdf.(Normal.(X[1: N-1], σ_x), X[2:N]))
# 観測データに関する対数密度
@views observation(X, Y, Z, w, σ_y) =
sum(logpdf.(Normal.(w[1].*Z .+ w[2] + X, σ_y), Y))
# 対数同時分布
log_joint_tmp(X, w, Y, Z, σ1, σ_x, σ_y, σ_w, N) =
transition(X, σ1, σ_x, N) +
observation(X, Y, Z, w, σ_y) +
prior(w, σ_w)
@views log_joint(X_vec, Y, Z, σ1, σ_x, σ_y, σ_w, N) =
log_joint_tmp(X_vec[1:N], X_vec[N+1:N+2], Y, Z, σ1, σ_x, σ_y, σ_w, N)
# log_joint(X_vec, Y, Z, σ1, σ_x, σ_y, σ_w, N) =
# transition(X_vec[1:N], σ1, σ_x, N) +
# observation(X_vec[1:N], Y, Z, X_vec[N+1:N+2]) +
# prior(X_vec[N+1,N+2], σ_w)
params = (Y_obs, Z_obs, σ1, σ_x, σ_y, σ_w, N)
# HMCの実行
x_init = randn(N+2)
maxiter = 1000
samples, num_accepted = inference_wrapper_HMC(log_joint, params, x_init, maxiter=maxiter, L=100, ϵ=1e-2)
println("acceptance rate $(num_accepted/maxiter)")
acceptance rate 0.964
# 推定結果の可視化
fig, axes = subplots(5, 1, figsize=(8,15))
# 出力値
axes[1].plot(Y_obs)
set_options(axes[1], "time", "y", "output data (Y)")
# 入力値
axes[2].plot(Z_obs)
set_options(axes[2], "time", "z", "input data (Z)")
# 回帰によって説明される成分
for i in 1:maxiter
w1, w2 = samples[N+1:N+2, i]
axes[3].plot(w1*Z_obs .+ w2, "g-", alpha=10/maxiter)
end
set_options(axes[3], "time", "w1x + w2", "regression")
# 状態変数によって説明される成分
for i in 1:maxiter
X = samples[1:N,i]
axes[4].plot(X, "g-", alpha=10/maxiter)
end
set_options(axes[4], "time", "x", "time series")
# 回帰と状態変数の和
for i in 1:maxiter
w1, w2 = samples[N+1:N+2, i]
X = samples[1:N, i]
axes[5].plot(w1*Z_obs .+ w2 + X, "g-", alpha=10/maxiter)
end
axes[5].plot(Y_obs, "k", label="Y_obs")
set_options(axes[5], "time", "w1x + w2 + x", "regression + time series", legend=true)
tight_layout()
